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ABSTRACT 

Detailed understanding of heat transfer and fluid 
flow is required for many aerospace thermal systems. 

These systems often include phase change and operate over 
a range of accelerations or effective gravitational 
fields . 

An approach to analyzing such systems is presented 
which requires the simultaneous solution of the 
conservation laws of energy, momentum, and mass, as well 
as an equation of state. The variable property form of 
the governing equations are developed in two-dimensional 
Cartesian coordinates for a Newtonian fluid. 

A numerical procedure for solving the governing 
equations is presented and implemented in a computer 
program. The Galerkin form of the finite element method 
is used to solve the spatial variation of the field 
variables, along with an implicit Crank— Nicolson time 
marching algorithm. Quadratic Lagrangian elements are 
used for the internal energy and the two components of 
velocity. Linear Lagrangian elements are used for the 


i 



pressure. 


The location of the solid/liquid interface as well as 
the temperatures are determined from the calculated 
internal energy and pressure. This approach is quite 
general in that it can describe he&t transfer without 
phase change, phase change with a sharp interface, and 
phase change without an interface. 

Analytical results from this model are compared to 
those of other researchers studying transient conduction, 
convection, and phase change and are found to be in good 
agreement. The numerical procedure presented requires 
significant computer resources, but this is not unusual 
when compared to similar studies by other researchers. 
Several methods are suggested to reduce the computational 
times. 
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NOMENCLATURE 

a thermal diffusivity 

C specific heat 

(also capacitance matrix) 

D geometric dimension 

e internal energy 

Fo Fourier number 

g acceleration 

(also gravity) 
i node number 

j time step increment 

Ja Jakob number 

k thermal conductivity 

K stiffness matrix 

L latent heat 

n unit normal vector 

m unit tangential vector 

M mass matrix 

N interpolation (shape) functions 

Nu Nusselt number 

n number of global nodes 

P pressure 

Pr Prantl number 

q heat flux 

r number of nodes per element 

Ra Rayleigh number 

T temperature 

t time 

u velocity in x-direction 

v velocity in y-direction 

V velocity 

W Gauss-Legendre weighting factors 

x coordinate in cartesian system 

y coordinate in cartesian system 

(i coefficient of thermal expansion 

Tj transformed coordinate 

p absolute viscosity 

p density 

a normal stress 

t shear stress 

(also dimensionless time) 
v kinematic viscosity 

0 field variable 

$ field value at a node 

C dimensionless location of phase interface 

4 transformed coordinate 
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parameter in time-marching recursion algorithm 
(also dimensionless temperature) 

¥ prescribed nodal value 

V divergence of a vector 


Subscripts and superscripts 

Q classifies an area or volume integral 
I classifies a surface integral 

* classifies an approximate value of field variable 

1 classifies liquid state 

s classifies solid state 

o classifies initial condition or reference state 

Matrix notation 

L J single row matrix 
{ } single column matrix 
[ ] matrix 
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CHAPTER 1 


INTRODUCTION 


The need for general analytical tools for modeling 
heat transfer and fluid flow is increasing as man designs 
more complex thermal/fluid devices. This is particularly 
true in the aerospace industry where highly reliable 
systems must operate in environments where little or no 
supporting experimental data is available. Such systems 
often include phase change and operate, over a range of 
accelerations or effective gravitational fields. 
Experimental investigations of fluid/thermal systems under 
low gravity conditions are difficult and expensive. 

Because of the time required for many phase change 
problems, most experimental studies are not possible in 
ground-based low— gravity facilities and must be done on 
Earth-orbiting laboratories. For these reasons a 
predictive analytical or numerical method would be very 
valuable. 

The intent of this research is to develop a general 
purpose numerical approach and computer program for 
analyzing the heat transfer and fluid flow of materials 
undergoing phase change. Such an analytical tool would 
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significantly reduce the number of experiments required 
and aid in our understanding of the experimental results. 

Many practical applications of such a computational 
tool exist, such as modeling cryogenic fluid management 
systems and analyzing advanced material processing and 
casting methods. Systems such as batteries and thermal 
management devices could be examined as they might have to 
withstand initial or inadvertent freezing in the low 
temperatures of space. Another application is in the 
analysis of designs for a thermal storage device to be 
used in the space power system on the National Aeronautics 
and Space Administration (NASA) Space Station Freedom or 
lunar base. As described by Klann 34 , this space power 
system would collect and concentrate solar energy to heat 
the working fluid of a Brayton cycle heat engine. A 
latent heat thermal storage device would provide energy 
for the power system during dark periods of the Space 
Station orbit, (see Burns 10 ) . 

Analytical approaches and numerical techniques for 
modeling thermal and fluid problems have been the subject 
of research and development for many years. Today, 
computer programs for modeling heat transfer by conduction 
are well developed and generally easy to apply. Until 
recently, computer programs for modeling fluid flow and 
its effect on convective heat transfer, however, had 
mostly been limited to empirical relationships based on 
known and simple geometries. Today, computational 
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convection using numerical methods such as finite 
differences and finite elements have become available to 
handle more complex flow geometries. Some of these are 
commercially available, however most are research oriented 
and limited in scope to particular applications. The 
additional complexities of having phase change phenomenon 
and materials that exhibit widely varying properties 
restricts the application of most present methods for such 
problems. Hence further work is needed to develop general 
purpose methods to analyze fluid/ thermal problems with 
phase change. 

In Chapter 2 a review is presented of the literature 
dealing with numerical solutions to thermal, fluid flow, 
and phase change problems. Chapters 3 and 4 cover the 
development of the governing equations and the numerical 
approach. Results and verification of the approach and 
model are presented in Chapter 5. The main body of the 
thesis ends with Chapter 6, in which concluding remarks 
and recommendations for further work are given. A 
computer program, PHASTRAN, developed during this research 
is discussed and presented in the appendices. 



CHAPTER 2 


LITERATURE REVIEW 


Much research has been devoted to the analysis of 
materials undergoing phase change because of its 
association with many applications. The food, 
metallurgical, and semi-conductor industries are important 
examples. More recently, there has been an interest in 
modeling these processes in the space environment. 

Many examples of research are in the literature for 
modeling phase change and fluid flow. For a detailed 
discussion on related subjects the reader is referred to: 
Stefan and Lunardini , on the phase change problem; 
Carnahan, Luther, and Wilkes 11 on numerical methods 
including the finite difference method; Baker 4 , Huebner 
and Thornton 12 , Zienkiewicz^ on the finite element 
method; Arpaci and Larsen 2 on convective heat transfer; 
and VanWylen and Sonntag^ 1 on thermodynamics. 

To summarize the some of the most important works 
related to this research topic, first those papers related 
to numerical modeling of the phase change problem will be 
discussed followed by those related to the fluid flow. 

The classical analytical approach to the phase change 
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problem will not be discussed. Though elegant in its 
mathematical derivation, its application to problems of 
complex geometry and temperature dependent material 
properties is impractical . Some noteworthy papers on the 
classical analytical modeling of phase change include 
Budhia and Kreith 9 , Siegel and Savino 54 . 

Most of the numerical models of phase change have 
involved the finite difference method, and to a lesser 
extent, the finite element method. This is not because of 
any superiority of the finite difference method, but 
rather the chronological development of the two methods. 
Indeed, most of the recent works have been devoted to the 
finite element approach. 

Otis 46 solved the melting problem by dividing the 
region into finite space intervals. Temperature was 
assumed uniform within each volume element at any instant 
and the latent heat effect was modeled as a moving heat 
source. The method required a coordinate transformation 
in terms of a pseudo time variable and was limited to 
analyses of materials initially or finally at the fusion 
temperature . 

Murray and Landis 4 ^ - suggested an approach by which 
the interface location was calculated by solving a 
differential eguation for the velocity of the phase front. 
The differential equation was derived from an energy 
balance at the phase front. The temperature at the new 
front location was then set equal to the freezing 
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temperature . 

Springer and Olson 56 used the Murray and Landis 
approach to track the phase front in two dimensions. 

Again the temperature at the phase front was set equal to 
the fusion temperature and the temperatures in the 
remainder of the solid and fluid was determined from a 
finite difference solution for heat conduction. 

Shamsundar and Sparrow 53 used enthalpy as the 
dependent variable instead of temperature in a finite 
difference formulation. Because their formulation 
involved an integral approach to the energy balance, the 
method eliminated the need to explicitly track the 
interface. They maintained this was the best method for 
analysis of multidimensional conduction phase change. 

More discussion on this method follows in Chapter 3. 

Only a few researchers have included the effects of 
natural convection in the fluid or radiative heat 
transfer. Such effects introduce nonlinearities in the 
field equations which require iterative solution 
procedures and increased computational times. 

, RQ 

Tien 3 solved the phase change problem with natural 
convection included in the fluid. He used a finite 
difference formulation of the conservation laws using a 
vorticity and stream function form of the momentum 
equations. Again the Murray and Landis approach was used 
to track the phase front. Tien's numerical results 
compared favorably with experimental data on the freezing 
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of naphthalene. 

Valle 60 also included natural convection in his 
solution but solved the problem using the finite element 
method. The conservation laws were developed in terms of 
the stream function and temperature. The latent heat 
effects and phase front motion were formulated implicitly 
in terms of an imbalance of the heat fluxes at the 
solid/liquid interface. This was one of the most detailed 
analyses of the phase change problem to date and included 
fluid flow, surface tension, and radiation effects. Valle 
compared his results to the work of Tien, however, and 
concluded that this approach did not seem to track the 
interface motion as effectively as approaches based on 
that of Murray and Landis. 

More recently, several works have approached the 
phase change problem using moving and deforming finite 
element grids and/or coordinate transformation. 

Ettouney and Brown 16 transformed the problem so that 
the melt and solid regions have fixed boundaries, of which 
the interface is one. This is an elegent approach which 
couples the interface shape and field variables allowing 
more efficient solution techniques. However, this 
approach, as with other moving mesh formulations, has the 
limitation of not being able to easily handle 
disappearance, merging or fragmentary distribution of 
phases . 

Albert and O'Neill 1 used a method of transfinite 
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mappings in conjunction with a moving boundary-moving mesh 
finite element technique. Improvement in tracking the 
phase front was made compared to a fixed mesh approach. 
Again there is the limitation mentioned above for the 
Ettouney and Brown method which restricts its application. 

Because of the high computational costs associated 
with modeling the phase change problem, some researchers 
have studied less numerically intensive schemes. 

Schneider formulated the phase change problem using the 
finite difference technique along with a variation of the 
enthalpy method of Shamsundar and Sparrow^"* ’ Depending on 
the amount of movement of the interface, Schneider's 
algorithm adjusts the number of convergence iterations. 

If the interface only moves within one grid spacing, only 
one iteration is used to converge the nonlinearities. 

This significantly reduces the computational times but may 
also affect the accuracy, especially for materials with 
properties that vary rapidly near the interface. 

The application of numerical methods to the modeling 
of fluid flow problems has made remarkable progress over 
the last 25 years. Initially, computer-based solutions 
used the finite difference method. Over the years, the 
finite difference method has provided solutions to many 
difficult flow problems including slow viscous flows, 
boundary layer flows and even variable property flows 
(thermo-hydrodynamic) flows. More recently, the finite 
element method has been developed to handle many of the 
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same problems. The finite element method has been shown 
to be particularly useful in problems involving complex 
geometries and boundary conditions. Baker , and 
Gallagher, et al . 21 contain many examples of the 
application of finite element to complex problems. 

Early applications of the finite element method to 
some continuum problems often used variational methods to 
derive the finite element equations. The lack of exact 
variational forms of the Navier-Stokes equations, however, 
prevented the use of finite elements to practical flow 
problems. Later, the application of weighted residual 
methods broadened the application of finite elements to a 
variety of fluid problems. 

Olson^ applied a pseudo-variational approach to a 
two-dimensional incompressible formulation developed in 
terms of the stream function. 

Baker 2 applied the weighted residual technique of 
Galerkin 22 to viscous incompressible flow. The Galerkin 
criteria, originally a nondiscretized approach is 
currently the most widely used method of formulating the 
finite element (discretized) equations. 

Hood and Taylor 31 also used the Galerkin criteria and 
formulated the Navier-Stokes equations in three ways: the 
velocity/pressure formulation; the stream function and 
vorticity formulation; and the purely stream function 
formulation. Comparison of these three formulations 
suggests that the velocity/pressure formulation may have 
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several advantages. It is readily extended to three 
dimensions. Pressure, velocity, velocity gradient, and 
stress boundary conditions can be easily handled. And it 
appears to reguire less computational time than the other 
formulations. 

Recently, more attention has been given to the 
considerations for obtaining good quality solutions to 
these nonlinear fluids problems over a wider range of 
conditions. Important aspects of this include proper 
choice of solution technique, element types, and mesh 
refinement. 

Gartling, et al . 23 ' formulated the finite element 
equations in terms of velocity and pressure and studied 
the convergence properties of several solution algorithms, 
two element types, and several mesh refinements. Laminar 
flow between converging plane walls was used to represent 
a nonlinear problem. Of the solution techniques, they 
found that those which solved the full unsymmetric 
equation system were superior and more generally 
applicable than their symmetric counterparts. In 
particular, the Newton-Raphson procedure was the most 
rapidly convergent. No significant difference was found 
between an 8-node quadrilateral element and a 13-node 
quadrilateral element, with the 8-node being preferred 
because of its reduced complexity in formulation and use. 
Finally, adequate mesh refinement was required in the 
direction of most rapid variation of the solution field. 
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Ben-Sabar and Caswell 6 investigated the effect of the 
choice of boundary conditions on the problems where the 
ratio of convective to diffusive terms are large. They 
found that consistent use of the velocity and surface 
traction boundary conditions are necessary to delay the 
appearance of numerical instabilities with increasing 
Reynolds number. 

Fletcher 20 developed an alternating direction 
implicit finite element method for flows where the 
convection terms dominate and applied the method to 
viscous compressible flow past a rectangular object. In 
comparison with an equivalent finite difference scheme, he 
found the finite element approach to be computationally 
more efficient. 

Solutions to coupled fluid/thermal problems continue 
to be the subject of much research, particularly transient 
problems in three-dimensional space. Such problems often 
require the expenditure of significant computer resources. 
Though beyond the scope of this study, a number of efforts 
are directed at improved solution methods to solving large 
systems of nonlinear equations. The use of many of these 
methods will be dependent on the availability of new 
computer architectures providing vector and parallel 
processing capabilities. 


CHAPTER 3 


PROBLEM FORMULATION 
Problem Statement 

The problem selected is to analytically determine the 
transient temperatures, heat transfer rates, fluid 
velocities, and pressures in a pure substance or eutectic 
material undergoing phase change. The material can exist 
in solid and fluid states with variable properties 
satisfying a general equation of state model. It is 
contained in a vessel of arbitrary geometry such as is 
shown in Figure 3.1. Boundary conditions could include 
prescribed temperatures, heat flux, and fluid velocities. 
Flow in the fluid is induced due to a gravitational body 
force or accelerating reference frame and both inertial 
and viscous effects are included. 

No surface free energy (surface tension) effects are 
included, nor is heat transfer by radiation. The effect 
of supercooling and mechanisms of nucleation or 
crystallization are also not considered. In addition, the 
fluid motion is restricted to laminar flow. 
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Figure 3.1 Example Phase Change System 


Governing Equations 

Problems in science and engineering can be classified 
as Lagrang ian or Euler ian depending on the viewpoint or 
reference frame adopted. To formulate the governing 
equations, one of these two approaches must be adopted. 

In the Lagrangian approach, all matter consist of 
particles which can be identified as they move through 
space. The independent variables in the Lagrangian system 
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are x Q , y 0 , z 0 and t where x Q , y 0 , z 0 are the coordinates 
which a specified fluid element passed through at time t Q . 

In the Eulerian approach, processes are characterized 
by continua of field quantities. The independent 
variables are the spatial coordinates x, y, z and time. 

To derive the governing equations, we focus our attention 
on one area in space called a control volume. If we apply 
the governing laws of the problem to a differential 
control volume we obtain a set of governing differential 
equations. This is the approach with which most problems 
in fluid and thermal analysis are formulated and is the 
approach adopted here. 

The solution to modeling the phase change problem 
includes solving the equations expressing the three 
physical laws of: 

1. Conservation of energy 

2 . Conservation of momentum 

3. Conservation of mass 

as well as a thermodynamic equation of state. 

Because this problem is dominated by thermal aspects, 
it is particularly important to consider the form of the 
energy equation. The conservation of energy is most 
commonly expressed in terms of temperature and specific 
heat. Such formulations are quite valid for single phase 
problems, however, they may be inappropriate for materials 
undergoing phase change at a discrete temperatures. 

To further discuss this, let us consider two 
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situations, one in which a sharp interface is formed and 
the other in which a mushy region with no sharp interface 
will exist. 

An example of a sharp interface might be a thin layer 
of water with its top and bottom sides insulated as shown 
in Figure 3.2. Suppose the water was initially at a 


Initial Temperature, 




Figure 3.2 System Exhibiting a Sharp Interface 
During Phase Change 


temperature above the freezing point and then one end of 
it is reduced to a temperature below the freezing point. 
Under these circumstances, a sharp interface will form 
which will separate the solid and liquid regions. In a 
volume element containing the interface however, the 
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specific heat is not easily defined. For such multi-phase 
problems, the energy equation is usually written 
separately for the solid and liquid phase regions. Since 
the interface is generally of unknown shape and position, 
numerical methods such as finite difference or finite 
element discretization encounter significant problems in 
handling the interface. Some numerical methods track the 
location of the interface. A heat balance can be 
formulated at the interface for more than one spatial 
coordinate. A differential equation is used to relate 
velocity of the interface to the heat absorption or 
removal. For three-dimensional phase change systems the 
interface is a surface and numerical methods for tracking 
the interface can become quite complicated. 

The temperature and specific heat formulation also 
has difficulty when analyzing phase change where the 
interface is not sharp. White® ^ justifies the existence 
of mushy regions with an example of the welding of two 
plates. A similar example would be to consider a 
thermally insulated plate with electrical connections at 
each end, as in Figure 3.3. When an electric current is 
passed through the plate it will heat up due to heat 
generation from internal resistance. Eventually, the 
temperature of the plate will reach the fusion (i.e. 
melting) temperature of the material. At this time the 
internal energy will equal the saturated value at that 
pressure and temperature and is equal to the product of 
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Electric 

Current 


Insulated Plate 



Uniform Heat Generation 
From Internal Resistance 


Figure 3 . 3 System Which Does Hot Exhibit a Sharp 
Interface During Phase Change 


the density, p, specific heat of the solid, C s , and fusion 
temperature, T f . This is shown graphically in Figure 3.4 
as point 1. As time continues, the internal heat 
generation results in increased internal energy with no 
change in temperature. This continues until the amount of 
energy converted to heat equals the latent heat and the 
material becomes liquid. This is seen as point 2 in 
Figure 3.4. As time progresses further, temperature of 
the material resumes its increase governed by the specific 
heat in the liquid. From this example, it is apparent 
that with the temperature description of the problem, the 
continuous transition is lost during the phase change. 
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Also assuming negligible diffusion of the thermal energy 
out of the electrodes, the material can exist at the 
fusion temperature in a two-phase state with no apparent 
sharp interface. 



Temperature 


Figure 3.4 Internal Energy Versus Temperature for a 

Substance which Changes Phase at a Discrete 
Temperature 


Many materials exhibit the behavior of phase change 
at a discrete temperature. To avoid the noncontinuous 
behavior of the product of temperature and specific heat, 
many formulations assume that the phase change occurs over 
a small but finite temperature range, see for example 
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Bamberger, et al . 5 This approach essentially defines an 
artificial specific heat for the volume containing the 
phase front. While these formulations retain temperature 
as the primary unknown, they may introduce significant 
errors in the results. Bonacina, et al.® demonstrated for 
example, that even in the one-dimensional case, the 
magnitude of the assumed range of phase change 
temperatures can affect the results significantly. 

After studying the above examples, it is apparent 
that it is the energy in a given volume that is really of 
interest and not its temperature. Thus, an alternative to 
the temperature and specific heat formulation is to use 
internal energy as the primary unknown and compute the 
temperature from the internal energy. - Shamsundar and 
Sparrow 53 were among the first to employ such an approach. 
They used an integral relation setting the rate of 
increase of the energy content in a arbitrary control 
volume equal to the net rate at which heat is conducted in 
through its surface. This relation was applicable whether 
or not the interfacial surface passes through the control 
volume. By assuming no fluid motion, pressure is 
independent of time, and they reformulated the problem in 
terms of enthalpy. Such a formulation turns out to be 
quite general in that it can describe heat conduction 
without phase change, phase change with a sharp interface, 
and phase change without an interface. Because the phase 
front is not explicitly tracked, the enthalpy formulation 
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avoids many of the numerical difficulties associated with 
fixed grid numerical methods, particularly in problems 
involving fragmented phases. 

Using an Eulerian frame of reference the governing 
equations for this problem can now be presented. For a 
variable property, Newtonian fluid, neglecting internal 
heat generation, surface free energy, and radiation, the 
conservation laws in Cartesian coordinates are: 


Conservation of energy : 




de 

+ 




+ f*r(x,y,t) 


(3.1) 


in which f is the viscous dissipation function given by 



For the problems presented in Chapter 5, the natural 
convective velocities are slow and the vicosities are low. 
Under these conditions, the viscous dissipation terms may 
be neglected. 
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Conservation of momentum (Navier-Stokes equations) : 



dU , 

pu a^ + 


3u 

pv a7 


pg_ 


ap 

ax 




(3.2) 
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+ pu ax + pv a7 


pg. 


ap 

ay 


do dT 



(3.3) 


in which 


o 
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2 p fx " f pVV 


(3.4) 


0 

y 


- f pVV 


(3.5) 


T = T 
x y y x 




(3.6) 


Conservation of mass : 


3 £ djyv _ 

at ax ay 


(3.7) 


Additional equations are required to evaluate the 
thermodynamic of state and material properties: 
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temperature, T=T(e,P) 
density, p=p(e,P) 
thermal conductivity, k=k(e,P) 
viscosity, p=p(e,P) 

If we can evaluate the state and the material 
properties explicitly we can reduce the number of 
equations to four and solve in terms of the basic unknowns 
u, v, P, and e. Specifying the state of a pure substance 
requires a minimum of two independent properties. When 
two phases of a pure substance exist together in 
equilibrium, the pressure and temperature are not 
independent and can therefore not be used to define the 
state. The two independent properties chosen as the basic 
field variables in the above equations are internal energy 
and pressure. 

The initial conditions consist of specifying the 
velocities, pressure, and energy at time zero. The 
hydrodynamic boundary conditions specify either the 
velocity components or surface tractions. The thermal 
part of the problem requires the heat flux or internal 
energy be specified on the boundary. Temperature boundary 
conditions must therefore be converted to internal energy. 

It should be noted that momentum and continuity 
equations as well as the convective transport terms in the 
energy equation are not required in that part of the 
solution domain which is in the solid state. The approach 
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to handling this problem is presented in the next chapter 
concerning the numerical method. 

At this point, many formulations, for example 
Valle**®, perform a transformation with the fluid velocity 
variables into a streamf unction and vorticity formulation. 
This was not chosen here. This decision was due to the 
requirement that this formulation be easily extendable to 
three-dimensional space. The streamfunction-vorticity 
formulation is often applied to two-dimensional 
incompressible flows. It can, however, be applied to a 
broader class of problems. This is because the 
definitions for the dependent variable transformations are 
essentially vector identities. These transformations can 
therefore be applied to three-dimensional Cartesian 
coordinate system. Unfortunately, for three dimensions, 
six scalar components for the streamfunction and vorticity 
must be defined compared to the four (3 velocities and 1 
pressure) used in the physical variable formulation. In 
addition, certain boundary conditions become difficult to 
apply. These difficulties have generally precluded 
application of streamfunction-vorticity to three- 
dimensional problems. 

As given in the problem statement, the present 
analysis is restricted to laminar flow conditions. Fluid 
motion is characterized as laminar if the fluid flows in 
imaginary layers and there is no macroscopic mixing 
between adjacent fluid layers. A flow is said to be 
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turbulent , however, if such nixing occurs. It should be 
noted that the governing equations given above hold at any 
instant of tine and apply to both laninar and turbulent 
flows. In a turbulent flow, however, the fluid velocities 
are fluctuating randonly about theif nean values. Such a 
randon variation in the field variables is nearly 
inpossible to solve directly. The standard approach is to 
tine-average the equations to obtain new ones which 
describe the tenporally averaged field variables. Such an 
approach is beyond the scope of the present formulation, 
and the flow is assumed laninar. 



CHAPTER 4 


NUMERICAL APPROACH 

The governing equations for the mass, momentum, and 
energy conservation given in the preceding chapter are 
represented by a system of nonlinear partial differential 
equations. These equations can describe some of the most 
interesting phenomenon in the fluid and thermal sciences. 
Unfortunately they are also some of the most difficult to 
solve. 

With few exceptions, (see for example Graebel ), 
problems involving convection can not be solved by direct 
integration of the partial differential equations. For 
most problems we must resort to some numerical solution 
method. In the approach used here, the finite element 
method is used to solve the spatial problem along with a 
recursive time marching algorithm based on the finite 
difference method. 

The finite element method is relatively new with most 
of its development occurring after 1960. There are 
several approaches to developing the finite element 
equations including the variational method, the method of 
weighted residuals, and the energy balance method. The 
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classical variational method is quite general, however it 
does require the existence of an exact variational form 
for the governing equations. For many problems, 
particularly in convective heat transfer, there are no 
exact variational forms. This requirement has limited the 
application of the variational method. Another procedure, 
the method of weighted residuals, or Galerkin method, does 
not require an alternate formulation of the physical 
problem and in fact can be applied to almost any well- 
posed system of differential equations. Oden 42 introduced 
a method by which the finite element equations can be 
developed from global energy considerations. This method 
has also proven very useful in the solution of many 
thermomechanical problems. 

Of the methods mentioned above, the Galerkin method 
has proven to be the most general and is the method chosen 
for the formulation developed here. 

The basic approach of the finite element method is to 
divide the solution domain up into a finite number of 
subdomains called elements. These elements are connected 
at node points on the element boundaries. The behaviors 
of the unknown field variables are then approximated 
within each element by continuous functions expressed in 
terms of nodal values of the field variables and their 
derivatives. Substitution of this approximation into the 
original differential equations and then integrating, 
results in some error or residual. In the Galerkin 
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method, linearly independent weighting functions are 
chosen such that the residual is required to vanish in 
some averaged sense over the entire solution domain. The 
resulting equations for each element are assembled into a 
set of coupled equations. 

The coupled equations are then directly integrated in 
time to yield the nodal values of the field variables. 

This direct integration of the coupled equations uses a 
recursion technique based on the finite difference method. 
Approximations in the finite difference method, however, 
introduce numerical errors. Though these errors can be 
minimized as the time step used approaches zero, it is at 
the expense of increasing computational time. Large time 
steps, in contrast, can produce entirely unrealistic 
behavior, including nonphysical oscillations which can 
even become unstable. Development of the proper recursion 
technique is thus of great practical importance and is 
discussed more, later in this chapter. 


Subd i v idina the Domain 


The first step in applying the finite element method 
is to subdivide the solution domain into elements. The 
selection of proper element type is still somewhat of an 
art. Lower order polynomial elements are simplest to 
formulate, but more elements are required for good 
solution accuracy. Fewer higher order elements are needed 
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for the same accuracy but require increased computation 
time in the numerical integration of each element. In 
general, to model a complicated boundary, it is usually 
more efficient to use a large number of simple elements 
rather than a few complex ones. Thus, for most problems, 
elements with interpolating functions of order greater 
than 3 are seldom used. 

In addition to computational efficiency, it is 
important that we select element types with interpolation 
functions that satisfy certain continuity and convergence 
requirements. This is necessary to ensure accuracy during 
integration and also that the approximate solution will 
converge to the correct solution with increasingly finer 
subdivisions (smaller elements) . These requirements were 
given by Felippa and Clough 19 and verified by Oliveira 43 . 
Specifically, they can be stated as 

1. The field variable <f> and its derivatives up to 
one order less than the highest-order derivative 
of the element (weak form) equations must be 
continuous at the element interfaces . 

2. The field variable 0 and its derivatives up to 
the order of the highest-order derivative of the 
element (weak form) equations must be continuous 
within the element . 

The first of these is known as the compatibility 
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requirement and the second as the completeness 
requirement. Compatibility requires that the field 
variable and its principal derivatives be the same at 
coincident nodes of neighboring elements. This ensures 
that there will be no contribution to the finite element 
equations from "gaps" at the element interfaces since the 
boundary integrals of each element will cancel. The 
completeness requirement ensures convergence to the 
correct solution when, in the limit, the element size 
shrinks to zero. 

It is convenient to introduce a standard notation to 
describe the degree of continuity of a field variable at 
the element interfaces. If the field variable is 
continuous at the element interfaces, it is said to have 
C° continuity. If, in addition, the second derivatives 
are also continuous, there is C 1 continuity, and so on. 

By choosing the internal energy, pressure and velocity 
form of the governing equations only first derivatives of 
the field variables appear. Thus only elements which 
satisfy C° continuity are needed to satisfy the above 
requirements . 

But other considerations may also influence the 
selection of proper element types. Several researchers 
modeling fluid flow have established that the 
interpolation functions for the velocity components must 
be at least one order higher than the pressure 
interpolation functions to prevent oscillations of the 
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field variable solution. Yamada, et al. 64 came to this 
conclusion by using a variational formulation. Olson and 
Tuann 45 showed that spurious rigid body modes in the 
solution appear when this criteria is violated. Other 
researchers who have supported this conclusion include 
Hood and Taylor 31 and Bercovier and Pironneau 7 . The 
restriction on the interpolation functions for the 
primitive variables arises from the uncoupled nature of 
the Navier-Stokes and continuity equations. This is 
because the continuity equation is simply a constraint on 
the velocities rather than an equation which fully couples 
velocities and pressure as the momentum equations do. A 
number of researchers modeling fluid flow using finite 
elements have concluded that quadratic interpolation 
functions for velocity and linear interpolation functions 
for pressure generally give the best performance 6 . An 
alternative approach to avoiding this problem is to 
uncouple the velocities and pressure by using a segregated 
method of solution. This is commonly done in finite 
difference formulations of the fluid equations, but it 
requires an additional convergence iteration to 
alternately satisfy the continuity and momentum equations. 

Two useful sets of rectangular elements are the 
serendipity and Lagrangian families. The serendipity 
elements shown in Figure 4 . 1 contain only boundary nodes 
and their interpolation functions were derived by 
inspection. The Lagrangian elements shown in Figure 4.2 




Bilinear Biquadratic Bicubic 


Figure 4.1 The Serendipity Family of Elements 



Bilinear Biquadratic Bicubic 

Figure 4.2 The Lagrangian Family of Elements 

contain interior nodes and use the Lagrange polynomial as 
its interpolating function. Both the serendipity and 
Lagrangian element types have seen wide use in finite 
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element analysis. 

As mentioned earlier, the geometry of the problem may 
also influence the selection of the element. One approach 
to modeling complex arbitrarily shaped boundaries is to 
use a body-fitted coordinate system. This approach can 
however add significantly to the modeling complexity. An 
alternative is to use curve-sided elements. Isoparametric 
elements are particularly useful as curve-sided elements. 
Isoparametric elements are elements whose geometry and 
field variable representations are described by 
polynomials of the same order. Using curve-sided 
elements, significantly fewer elements are usually 
required to fit a complex geometric boundary. Curved- 
sided isoparametric elements are commonly formed from 
either serindipity or Lagrangian elements. 

Finally, other numerical considerations may also 
influence the selection of the proper element type. The 
numerical solution approach for the problem in this 
research requires the use of element mass lumping to 
prevent unrealistic oscillations in the field variables. 
This will be further discussed in the section on solving 
for the transient response. The use of element mass 
lumping has been shown by Gresho, et al . 28 to yield 
unstable solutions with the quadratic serindipity element 
under certain conditions. The Lagrangian biquadratic 
element however showed good accuracy and stability. 

During the course of this research, stability problems 
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were also observed when using the quadratic serindipity 
elements with element mass lumping. 

For the reasons given above, the 4-node Lagrangian 
linear element was chosen for the pressure field and the 
9-node Lagrangian biquadratic element was chosen for the 
energy and velocity fields. 

Besides the element selection, the subdivision of the 
domain can have a significant influence on the solution. 

It is easiest to generate a uniform element mesh, however, 
this may not always provide the best representation of the 
field. Usually more elements should be placed in regions 
where the boundary is irregular. Also, in general, the 
elements used should be well proportioned, with the ratio 
of their largest dimensions to their smallest dimensions 
near unity. Nevertheless, it can be acceptable to use 
long thin elements if it is known that the field does not 
vary greatly in the elements lengthwise direction. 

Provided that elements have been selected which 
satisfy the compatibility and completeness requirements, 
increasing the number of elements will provide improved 
solution accuracy. If there is an approximate solution to 
the problem, the finite element model accuracy can be 
improved by using a finer mesh in areas where high 
gradients are expected in the field variable. This 
increased accuracy is at the obvious expense of increased 
computational effort. It is generally good practice to 
obtain several solutions to a problem using an increasing 
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number of elements. By comparing results it can then be 
determined what is a sufficient number of elements for 
good solution accuracy. 

Once the element type has been chosen, the 
interpolating functions for both the linear and quadratic 
elements can now be developed. The Lagrange polynomial is 
defined by 


L (x) 



(4.1) 


(x-x ) • •• (x-x ) (x-x ) 

o k - 1 _ x k + 1 ' 


(x-x ) 


(x -X ) 
k o' 


(X -X ) (X -X ) 
' k k - 1 ' ' k k + 1 ' 


(X-X ) 
k n 


Using the 4 -node rectangular element and local coordinate 
system defined by Figure 4.3 the variation of some field 
variable 0 can be written as 


<t>U. n) = N x + N 2 (4,n)* 2 




+ N (£,7))* 

4 4 


where # represents the nodal values of the field variable 
and the interpolating functions N are given by 


N x (4,r?) = L (4)L i (r,) , 


n 2 (4, r,) = l 2 (4)L 2 (n), 


etc. 
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Figure 4.3 The 4-node Lagrangian Element 
and Coordinate System 


These interpolation functions, formed as products of the 
Lagrange polynomial, are bilinear. The explicit 
expression for node 1 follows: 


N x U,n) = L (OL (n) 


{ - { 2 „ 


-1 

1-1 


-QZ 1 

-i-i 


In this manner all of the linear Lagrangian interpolation 
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functions can be developed and are given by 


N = J Uri-i-n+1) 
N 3 = J (^rj+4+n+i) 


n 2 = J (- 4 n+ 4 -n+i) 

N (^rj+4-rj-l) 


(4.2) 


The quadratic Lagrangian interpolation functions can be 
developed similarly and are given by 


n x - J Un~i 2 n~4n z +4 z n 2 ) 
N 3 = 4 (-£n- 4 2 n+ 4 r) 2 + 4 V) 
n 5 = "4 (iv+i 2 ti+4i 2 +4 2 n 2 ) 
n 7 = J (- 4 n+ 4 2 n- 4 n 2 + 4 2 n 2 ) 
n = i-$ 2 -n 2 +4V 


N 2 = | C-T?+T? 2 +4 2 n-4 2 1? 2 ) 

n = | ( 4 + 4 2 - 4 n 2 - 4 V) 

(4.3) 

N 6 = 2 

N 8 = 2 (“*+« 2 + * , » 2 -«V> 


After selecting the element type, the solution domain 
is subdivided into a specified number of these elements. 
Fitting a curved boundary such as shown in Figure 4.4 
could be done with many small elements, however in this 
case, a better fit would result if we could use curve- 
sided elements as in Figure 4.5. Ergatoudis et al. 17 were 
among the first to develop a general approach to creating 
such elements. Curved-sided elements are developed by 




Figure 4.5 Curve-sided Elesent 
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transforming or mapping simple geometric shapes in some 
local coordinate system ( 4 -rj) into distorted shapes in the 
global Cartesian coordinate system (x-y) . To construct a 
typical element such as is shown in Figure 4.5 we must 
start with a simpler "parent ' 1 element. Consider a parent 
element such as the 9 -node quadratic Lagrangian element 
shown in Figure 4.6. The coordinates in the plane may 
be transformed into the x-y plane using mapping functions 
of exactly the same form as the interpolation functions. 
These are given by 


9 9 

x= and y= ^^(4,1,^ (4.4) 

i =1 i =1 


When making this transformation we must of course 
ensure that for every point in the local 4-17 coordinate 
system there is a unique corresponding point in global x-y 
coordinate system. If the transformation is not unique, 
the element can be greatly distorted causing unpredictable 
results on the solution. 

This transformation technique can be useful in 
generating a set of element coordinates for a region in 
the solution domain such as that shown earlier. The 
region coordinates in the global cartesian system would 
become the x^ and y^ in equation (4.4). For a division of 
nine elements in the "parent" region shown in Figure 4.7, 



39 



Figure 4.6 The 9-node Quadratic Lagrangian Element 


the interpolation functions are evaluated at the 

appropriate 4 and n having discrete values of -1, -0.667, 
0.333, 0, -0.333, -0.667, and 1. The resulting elements 
in the global cartesian system are shown in Figure 4.8. 
This technique is quite useful for automatic element 
(grid) generation and is not restricted to equal numbers 
of divisions in the 4 and *1 directions. It is also 
possible to make the grid mesh finer in an area by slight 
shifts in the discrete 4 and t? values given above. For 
further discussion on grid generation techniques see 
Zienkiewicz, et al. 65 . 

Before substituting the interpolating functions into 
the finite element equations, it is also necessary to 
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develop expressions for their derivatives. Following 
Huebner 32 , the variation of some field variable <j> within 
an element having r nodes is again expressed as 


r 

<PU,ri) = ^ N (4,rj) 

i=l 


(4.5) 


The derivatives of the field variable can also be 
evaluated by 


ax 2— ax *i 

i=l 


and 



(4.6) 


To evaluate the element matrices we must also express 
3^/ax and 3^/ay in terms of local coordinates 4 and 77 . 
Applying the chain rule of differentiation yields 
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(4.7) 


where [J] is the Jacobian matrix. The Jacobian is 
evaluated using 




Figure 4.8 Curve— sided Eleients in the Cartesian Systen 
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[j(£,n)] = 


9N 9N 

/ H ( ^ ,r?)x i Z- df ^ ,r,)Y i 


i =1 

r 


i=l 

r 


r — 1 3N r- — • dN 

2_^‘(«.n)x 2 _ 5if l ({>i)yj 


i =1 


i =1 


( 4 . 8 ) 


Rearranging, the derivatives of the shape functions 
in the two coordinate system are related by the inverse of 

the Jacobian as follows: 
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r l " 1 
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l 
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for i = 1,2, .... ,r 


( 4 . 9 ) 


From the above equations we can find the partial 
derivatives of the field variable in terms of the 
transformed coordinates 4 arid rj using 
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( 4 . 10 ) 


Finally, from advanced calculus we can express the 
differential area dx dy in terms of d£ dr) using 



dx dy = |j| d^ drj 


(4.11) 


in which |j| is the determinant of the Jacobian. 

The validity of the element equations depends on the 
existence of the inverse of the Jacobian for each element. 
Also, the 4 -T J to x-y coordinate mapping discussed earlier 
is unique only if the inverse of the Jacobian exists. A 
useful method for determining this uniqueness and the 
validity of the mapping is to evaluate the determinant of 
the Jacobian for all elements. If the sign of the 
determinant does not change throughout the solution 
domain, an acceptable mapping will be assured. 

Formulating the Element Equations 

To formulate the finite element equations from the 
governing equations we must apply the Galerkin method, 
substitute the interpolation functions for the field 
variables and their derivatives, and then perform the 
numerical integration on an element basis. 

The velocity, pressure, and energy distribution 
within each element can be approximated by 


u(x,y,t) = Ln u ( x,y)J (u(t) ) 


(4.12) 
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v(x,y,t) = |_N v (x,y)J {v(t) } (4.13) 


P(x,y,t) = Ln p (x,y)J (P(t)} (4.14) 


e(x,y,t) = Ln* (x,y)J (e(t) ) 


(4.15) 


Before applying the Galerkin method to the 
conservation of energy equation it is necessary to first 
linearize the nonlinear convective terms. Let u* and v* 
be an approximate solution to the velocity field and P* be 
an approximate solution to the pressure field (such as the 
results from a previous iteration) . Now applying the 
Galerkin method, the linearized energy equation yields 



(4.16) 


Integrating the last two conduction terms by parts 


using Green's theorem 
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f a(V b) dn = f a(bn) dX - j b-Va dQ (4.17) 

J Q J I J n 


am ^ HT 

letting a = N* b = k a£ n i + k dy n j 


(4.18) 


the energy equation becomes 



The surface integral in the equation above is the natural 
boundary condition and allows for the introduction of the 
prescribed heat flux, q, boundary condition. 

Substituting the approximations for the field 
variables and rearranging, the energy equation for the 
interpolating function at node i becomes 
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( 4 . 20 ) 


Applying the Galerkin method now to the momentum 
equations it is necessary to also linearize the nonlinear 
convective terms, by again letting u* and v* be an 
approximate solution to the velocity field. Taking the x- 
direction momentum equation and applying Galerkin 's 
criterion yields 



* du 
pU 3 ^ 


* du 
pV 3 ^ 


3(o -P) 3 t 

x _ ij 
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+ pg 


( 4 . 21 ) 
dfl = 0 


Integrating the viscous force terms by parts using Green's 
theorem 


b = o n + r n 

x i xy j 


letting 


a = N v 

i 
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yields 


I n N I (a/' ’ + 37") dn - J/I (<•.-*>*» + T „ a ,) dI 

U 8N 9N ^ 


( 4 . 22 ) 


Now defining 


0 = (a -P) n +t n 

x ' x ' i x y j 


( 4 . 23 ) 


and introducing the velocity component's with 



fnv-v 



fjivv 


( 4 . 24 ) 
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The x-momentum equation becomes 
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( 4 . 25 ) 


Substituting the approximations for the field variables 
and rearranging, the x-momentum equation for node i of the 
element becomes 



( 4 . 26 ) 


The integral over the surface X is the natural boundary 
condition and can be used to introduce surface tractions. 
The y momentum equation can be developed similarly 


and is given by 
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( 4 . 27 ) 


Finally, applying the Galerkin method to the 
conservation of mass equation (continuity) , using 
weighting factors equal to the interpolating functions for 
pressure at each node i yields 



3£V] 

ay J 


dQ 


o 


( 4 . 28 ) 


Substituting the approximations for the field variables u 
and v yields 

J Q »; HJ *!«}♦/„< H5I d0 H - 0 < 4 - 29 

It is important to note that all three conservation 
equations are applied throughout the solution domain. 
However, since there is no motion of the material in the 
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solid something must be done to prevent convective 
transport and influence of mass on the buoyant forcing 
function. 

The simplest method for preventing fluid motion in 
the solid is to set the velocities u* and v* to zero and 
to increase the viscosity to some large value for those 
nodes which are in the solid state. As discussed earlier, 
during the phase change process, some materials develop a 
mushy region in which the interface between the phases is 
not sharp. If information is available on how the 
viscosity (and other properties) of the material varies as 
it undergoes phase change, it can be used in the material 
property data to improve the modeling of the flow near the 
phase front. For the cases studied in this research, the 
phase front was relatively sharp and the exact value of 
the assumed viscosity in the two phase region had little 
influence on the results. Also, because the results from 
the fluid equations are not applicable within the solid 
region, the exact influence of the viscosity in that 
region is unimportant. 

The second problem alluded to above involves the 
influence of the density distribution within the solid on 
the overall buoyant forcing function. In reality the body 
forces on the solid are balanced by internal and boundary 
stresses. In this formulation, however, no such mechanism 
exists since no equations from solid mechanics were 
included. Because of this problem, incorrect body force 
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terms can develop if the density of the solid is 
significantly different from that of the liquid. To 
prevent this, the buoyant force term is first reformulated 
in terms of density change about a reference density such 
as 



dn 


This approach is quite common in formulations of the 
Navier-Stokes equations which use the coefficient of 
thermal expansion. The reference density p 0 is taken to 
be the density of the fluid at the fusion temperature. To 
prevent the influence of the solid in the buoyant force 
terms, the density is set equal to the reference density 
for those nodes which are in the solid region. Note that 
this is done only in the evaluation of this integral while 
all other integrals are evaluated using the appropriate 
density for each state. 

Because the fluid motion is influenced by buoyancy 
forces and the material properties vary with time, the 
energy, momentum, and mass equations are directly coupled 
and must be solved simultaneously. Two approaches have 
been used in the past to solve the steady solutions. 

Taylor and Ijam®^ solved the equations simultaneously. 
Gartling 24 used an algorithm in which the equations are 
segregated and the solution alternates between the them. 
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During the course of this research, both approaches were 
employed and evaluated. The method which solves the three 
equations simultaneously was found to require 
significantly more computer memory and computations in 
solving the equations. The alternating solution method 
required an iterative algorithm, however, this did not 
substantially change the overall algorithm. This is 
because iterations are required to satisfy the nonlinear 
terms in the Navier-Stokes equations. Based on these 
studies, the segregated approach to solving the energy and 
flow equations was adopted in the present analysis. 

The energy, momentum, and mass equations above were 
given for the weighting functions at each node i in the 
element. By inspection, we can write -the energy equation 
for all nodes of each element as 


mi*} + m m 


R q + R p + R c J (4.30) 


Similarly, we can write the momentum and continuity 
equations for all nodes of each element as 
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To evaluate the finite element matrices requires 
integrating functions of the form 


f y) an = f f /'(4,rj) |j| d4 dr, 

J Q J -i J -i 


(4.32) 


The Jacobian is a function of 4 and r, and cannot be 
explicitly evaluated because the coefficients are 
polynomials, thus some type of numerical integration must 
be used. The Gauss-Legendre method is chosen here because 
it requires relatively few sampling points to obtain a 
good degree of accuracy. This method involves evaluating 
the function at the sampling points and weighting the 
results as follows 
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(4.33) 


Table 4.1 gives the location and weights for the Gauss- 
Legendre Quadrature up to order 4 
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Table 4.1 Location and Weights for Gauss-Legendre 
Quadrature to Order 4 


Order 

Location 

Weight 

n=l 

±0.5773502691 

1.0000000000 

n=2 

0.0000000000 

±0.7745966692 

0.8888888889 

0.5555555556 

n=3 

±0.3399810435 

±0.8611363115 

0.6521451548 

0.3478548451 

n=4 

0.0000000000 

±0.5384693101 

±0.9061798459 

0.5688888889 

0.4786286704 

0.2369268850 


Figure 4 . 9 shows an example of the location of these 
sampling points for a typical element using a Gauss- 
Legendre Quadrature of order 2. 

To accurately evaluate the volume integrals, Gauss- 
Legendre integration of order 1 is required for the 
bilinear elements and order 2 for the biquadratic 
elements. 


Assembling the System Equations 

Once the behavior of each element has been developed, 
the overall system is modeled by assembling these element 
equations into a set of system equations. To do this, the 
element equations, which were evaluated at the nodes of 
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Figure 4.9 location of the Gauss-Legendre Integration 
(Order 2) Sampling Points in a Typical 
Element 


each element, must now be transformed into the global node 
numbering scheme. The numbering schemes for the guadratic 
Lagrangian element and the global node numbering scheme 
for an example four element region are shown in Figure 
4.10. For the sake of explanation assume there is a 
single field variable at each node, the total number of 
system nodal variables is equal to the number of global 
nodes, n (i.e. 25 in Figure 4.10). The nine local element 
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nodes for the upper left element in Figure 4.10 correspond 
to the global nodes 11, 12, 13, 8, 3, 2, 1, 6, and 7. The 
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Figure 4.10 Transformation of the Local Node Numbering to 
the Global Node Numbering Scheme 


node relationship is only slightly more complicated for 
the case where multiple field variables exist at the same 
geometric node location. The procedure for assembling the 
system equations is as follows: 

1. For n global nodes, set up two n x n and one n x 
1 null matrices (all zero entries) as the system 
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mass, stiffness and resultant matrices. 

2. Take one element and use the relationship between 
the local and global node numbers to replace the 
indices in the element matrices with the 
corresponding global node numbers. 

3. Insert those terms into the appropriate locations 
in the system matrices. If a term is inserted in 
a location where another term has already been 
placed, it should be added to the value at that 
location. 

4. Repeat the procedure starting at step 2 for all 
of the elements. 

The result will be a system of equations of the form 


["]{♦} + [ k ]{♦} - !»<*>} 


(4.34) 


where again <f> is the unknown field variables, [M] the mass 
matrix, [K] the stiffness matrix, and {R} the resultant 
column matrix. 


Solving for the Transient Response 


The solution of the final set of simultaneous 
nonlinear ordinary differential equations is a formidable 
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task combining transient time integration with an 
iteration scheme at each time step to handle the nonlinear 
terms. 

An approach to directly integrating these coupled 
equations in time is to use recursive algorithms based on 
the finite difference method. Let t^ be a typical time in 
the transient response such that 

t =t +At , for j=0,l,2,... 

j + 1 j 

A general family of algorithms can be developed by 
introducing a parameter 0 such that 

t =t + 0At , for O^0^1 
i + e i 

The system equations at time t # can be written as 


[«](*u MR..- H.U 


(4.35) 


Introducing the approximations 
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Substituting these approximations into the system 
equations yields 


[HK] + [(«-!) [K] + £[«)]{♦}, 




(!-9)[r} j + •{«},„ 


(4.37) 


Rearranging, a general recursion formula for calculating 
the unknown field variables at the end of the time 

step to the known values {<f>) j at the start of the time 
step is given by 


[*]{*U- l s L. 


(4.38) 
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in which 

[ K ] - «[ k] + £[>] 

[( 9_i > i k i + AHHi + “->H* e H*i 

This equation represents a family of popular time-marching 
algorithms. Table 4.2 describes some of the members of 
this family. 


Table 4.2 Characteristics of Recursive 
Time-Marching Algorithms 


Algorithm 

e 

Accuracy 

Stability 

Euler or Forward 

0 

1st Order 

Conditional 

Difference 

Crank-Nicolson 

1/2 

2nd Order 

Uncond i t i ona 1 

Galerkin 

2/3 

1st Order 

Unconditional 

Backward Difference 

1 

1st Order 

Unconditional 


All of the algorithms given in Table 4.2 are first 
order accurate with the exception of the Crank— Nicolson 
method which is second order accurate. The terms first 
order and second order refer to the truncation errors in 
the finite difference approximations. First order 
accuracy means that the error is proportional to the first 
power of the time step At, and second order accuracy means 
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the error is proportional to the second power of At. 

Many studies have been made to evaluate the relative 
accuracy of each of the methods given in Table 4.2. 

Perhaps one of the most relevant was performed by Hogge 29 
in which he studied a nonlinear heat transfer problem. 

His detailed investigation of the relative accuracy of the 
various methods concluded that the Crank-Nicolson scheme 
(0=1/2) is indeed the most accurate of these methods. He 
did note however, that more sophisticated schemes spanning 
several time steps can give even better accuracy. 

Table 4.2 also characterizes the stability of the 
various methods. Stability means that the computed 
response does not oscillate and grow without bounds 
unrealistically. Stability is ensured- for- 0^1/ 2 . All of 
the methods given in Table 4.2 are unconditionally stable 
except the Euler forward difference method. For methods 
where 0<l/2 a stable solution results only for time steps 
less than some critical value. It should be noted however 
that the selection of time step is important even for 
methods with 0>l/2. Though the computed response with one 
of these methods will not grow unrealistically without 
bound, it may exhibit spurious oscillations and decreased 
accuracy with a very large time step. With either method 
it is good practice to solve the integration with several 
different time steps and compare the results. 

In addition to accuracy and stability, other 
considerations effect the selection of time integration 
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algorithm. The two common approaches to solving the 
overall system of equations are the explicit forward 
difference scheme and the implicit one parameter 0 
schemes . 

The explicit forward difference scheme computes the 
field variables at time tj+i from a set of uncoupled 
system equations. It does however require a lumped mass 
matrix. The term lumped is used to differentiate it from 
the original (or consistent) mass matrix. Lumped matrices 
are formed by assigning each node an amount of mass which 
can be attributed to that location. The most common 
approach to forming a lumped mass matrix is to sum the 
coefficients of the rows of the consistent mass matrix and 
use these as terms along the diagonal.- The explicit 
forward difference scheme using a lumped mass matrix may 
result in a significant computational savings over 
implicit schemes because the field variable can be 
computed without solving the system of simultaneous 
equations at each iteration. Again it has the 
disadvantage of only conditional stability with selection 
of the time step. It also requires a constant time step 
throughout the solution. 

The implicit M 0" algorithms compute the field 
variables from a coupled set of system equations. The 
time step for the implicit algorithms again is not 
restricted by the stability constraint discussed earlier. 
In addition, the time step can be varied throughout the 
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transient solution. This is a significant advantage of 
the implicit algorithms for this research problem. During 
the solution of a phase change problem, the maximum 
allowable time step may vary significantly. This is 
because the solution domain may be all solid, multi-phase, 
or fluid. Each of these situations could have very 
different rates of response. For example, during the time 
the material is all solid and only the energy equation is 
important, the allowable time step might be significantly 
larger than when fluid is present and results from the 
momentum equations become important. 

Implicit algorithms permit either lumped or 
consistent mass matrices. However, the choice of lumped 
versus consistent mass matrices is not always be easily 
resolved. Considering their formulation, consistent mass 
matrices are thought to be more accurate. Many 
researchers, however, have found insignificant loss of 
accuracy using the lumped approach. In fact, virtually 
a ll finite difference formulations use the lumped mass 
approach. Emery, et al . 16 found that the consistent 
approach could sometimes even predict unrealistic 
oscillations in the temperature distributions. This was 
most often observed near areas of sharp transients. The 
lumped approach however gave solutions which were 
intuitively obvious. During the course of this research, 

I also observed unrealistic oscillations in the field 
variables while using a consistent mass matrix approach. 
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This was particularly evident near the phase change front 
where the material properties varied greatly. 

Because of the considerations and observations 
discussed above, the lumped mass, implicit Crank - Nicolson 
scheme was adopted for this research. 

After application of the recursive time integration 
approach, the result will be a reduced set of system 
equations of the form 


MM -H 


It is now necessary to account for any boundary 
conditions which were not already applded as natural 
boundary conditions. In particular, these include any 
prescribed value boundary conditions. Usually, at least 
one and sometimes more than one nodal value must be 
prescribed to make the system equations nonsingular and 
provide a unique solution. There are several ways to 
apply these prescribed boundary conditions and modify the 
set of system equations. The one chosen here is 
relatively straightforward and can be described best by 
example. Suppose there are only four system equations 
given by 


68 


k k k k 

11 12 13 14 

k k k k 

21 22 23 24 

k k k k 

31 32 33 34 

k k k k 

41 42 43 44 



V 


r 

l 


*2 


r 

2 


*3 


r 

3 

. 

*4 


r 

4 


Consider applying prescribed nodal values specified as 
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The modified system equations will become 
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Note that once these prescribed boundary conditions 
have been applied, the number of equations and nodal 
unknowns to be solved for is reduced since it is not 
necessary to solve for the prescribed values. For 
problems with many prescribed boundary values, such as 
flow inside a containment vessel, the size of the system 
equations can be reduced significantly. From a 
computational standpoint, this can result in a substantial 
decrease in the time to invert the large system matrix. 

Note also that since the dependent variable in the 
energy equation is internal energy, prescribed temperature 
boundary conditions must be handled with an additional 
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step. This step is simply to apply the equation of state 
model to convert the prescribed temperature value to a 
prescribed internal energy value. 

Before summarizing the overall numerical approach, 
some discussion is warranted on the method for determining 
the state and material properties used throughout the 
solution domain. There were several requirements for this 
equation of state and material property model. First, the 
method should be able to describe reasonably complex 
material property characteristics. Second, the amount of 
input data to describe these material properties should be 
minimal. And finally, the method should be 
computationally efficient. Usually, either tabular data 
interpolation or "curve fitting" approaches are used for 
such models. Since two independent material properties 
are used (e.g. pressure and internal energy) the tabular 
data approach requires double interpolation and the 
"curve" in the second approach is really a surface. 

Because of the requirements discussed above, a 
surface fitting approach to the material state properties 
was developed. In such an approach, the thermodynamic 
surface representing the dependent property as a function 
of the two independent properties is represented by one or 
more regions as shown in Figure 4.11. These regions are 
described by quadrilaterals defined by values at eight 
points along their sides. The point representing the 
independent properties is projected onto one of these 
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Figure 4.11 Surface Fitting of the Material 
State Properties 


quadrilaterals. A double quadratic regression analysis of 
that quadrilateral is then used to yield the value of the 
dependent property. 

The overall numerical solution procedure described in 
this chapter is summarized in Table 4.3. This numerical 
solution procedure is implemented in the computer program, 
PHASTRAN, which is described in Appendix A and listed in 
Appendix B. 
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Table 4.3 Overall Solution Procedure 


Initial calculations 
Read input data 
Generate element coordinates 
Initialize field variables 
Evaluate material state properties 
At each time step 
Increment time 

Iterate to converge nonlinear terms 
Form the element equations 
Assemble [M] , [K] , and {R} system matric 
Modify system matrices to form [K] and { 
Apply the prescribed boundary conditions 
Solve the simultaneous equations 
Update the material state properties 
Set u* and v* to zero in the solid 
Check for convergence of the field variables 
Advance to the next time step 


Wlfl> 



CHAPTER 5 


RESULTS AND VERIFICATION 

Direct experimental verification of the multi- 
dimensional phase change problem is difficult at best. 
Verification is especially difficult for containment 
vessels with complicated shapes and nonuniform boundary 
conditions. Even for very simple geometries, verification 
would rely on results from other numerical methods or the 
very few experimental observations of t:he combined 
transient effects that do exist. For this research, an 
alternative but indirect approach was chosen in which the 
individual phenomena are verified with simple geometries 
for which there are known and well established solutions. 
This approach resulted in a number of test cases which are 
summarized in Table 5.1. These cases are presented 
individually in the remainder of this chapter. Table 5.1 
characterizes each case by the geometric space (1- 
dimensional or 2 -dimensional) of the problem, the 
thermodynamic state of the material, and the principle 
phenomena of interest . Note that even though the 
Particular case can be characterized as 1-dimensional, the 
2-dimensional analysis was used. 
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The last case gives the solution to the general 
problem combining all the phenomena of interest. This 
case demonstrates the ability of this analysis approach to 
solving a realistic problem representing several 
interacting fluid/thermal phenomena'. It could also serve 
as a test case for comparison of similar analyses which 
other researchers may be developing. 

A discussion of the computer resource usage is given 
at the end of the chapter. 


Table 5.1 Summary of Cases 


Case 

Space 

Material State 

Phenomena 

1 

1-D 

Solid 

Conduction, with 
Prescribed Temperatures 

2 

2-D 

Solid-Liquid 

Phase Change, by 
Conduction only with 
Prescribed Temperatures 

3 

2-D 

Liquid 

Buoyancy-Driven 
Convection with 
Prescribed Temperatures 

4 

2-D 

Solid-Liquid 

Phase Change , by 
Conduction and Buoyancy 
Driven Convection 
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gase 1; Conduction Only 

As was discussed in Chapter 3 , the formulation of the 
energy equation is based on internal energy instead of the 
more common temperature and specific heat approach. To 
ensure the validity of such an approach, Case 1 considers 
a simple transient problem for conduction heat transfer. 
For this case, material properties are constant throughout 
the solution domain and the material always remains in the 
solid state. Figure 5.1 describes the solution domain as 
well as the boundary and initial conditions. The problem 
consists of a slab of material initially at a uniform 
temperature, T o . At time zero the surface temperature of 
one side is suddenly changed to T t . Note that since the 
upper and lower sides are insulated, this problem is 
actually one dimensional. 

Exact solutions to the problem of Case 1 are widely 
available, for example Kreith 35 . The results of such 
solutions are often presented in terms of a nondimensional 
temperature versus the Fourier number defined by 

F = at 
° (2D) 2 

where a is the thermal diffusivity defined by 


75 


h«— 



Initial 

Temperature 

T=T n 


Figure 5.1 Description of Case l 



and D is the thickness of the material. 

Though this problem is really one-dimensional, for 
convenience, the solution domain was discretized into a 
total of 25 equal elements, with 5 element divisions along 
each of the x and y directions. The time step used 
corresponds to a Fourier number of 0.0005. Figure 5.2 
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Figure 5.2 Calculated Temperatures for Case 1 


shows the calculated temperature distributions at four 
times. These values are within one percent of the exact 
values for this problem. During the course of this 
research, many other similar problems involving conduction 
heat transfer were solved. These included cases with both 
prescribed temperatures and/or prescribed heat flux. 
Results from these cases as well as those of other 
investigators (see for example White® has confirmed the 
validity of the internal energy formulation of the energy 
equation in calculating transient conduction heat 
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transfer. 


Case 2: Phase Chai 


Conduction 


Several problems were studied involving phase change 
by conduction. For 1-dimensional space, exact solutions 
exist for prescribed temperature boundary conditions as 
well as cases with precribed heat flux conditions. 

Results from the present formulation for such cases showed 
excellent agreement with exact solutions. A more complex 
case is that of multi-dimensional phase change. Case 2 
models 2-dimensional phase change by conduction with 
prescribed temperature boundary conditions. Specifically 
the problem consists of a prism of square cross-section 
which is initially in the liquid state at the fusion 
temperature. Prescribed temperatures which are lower than 
the fusion temperature are applied to the surface of the 
prism, and it solidifies with time. Due to the symmetry 
of the problem, only one quarter of the cross-section need 
be considered with two boundaries maintained at the 
prescribed temperature as shown in Figure 5.3. For 
convenience the problem can be described by 
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Initially at Fusion 
Temperature, Tj 




Sudden Change 
in Surface 



Figure 5.3 Description of Case 2 


(VC ) 

k = t — t dimensionless latent heat parameter 

£ 1 


(T f - T ) 

0 f = "7^ — m > dimensionless fusion temperature 

' o 1 ' 


a t 

s 



T 


dimensionless time 
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Figure 5.4 Calculated Loci of Interface for Case 2 



dimensionless interface location 


in which D is a convenient reference length and the 
subscript s denotes properties of the solid. 

This problem has been studied by several researchers 
in the past including Poots 50 , Lazaridis 37 and Crowley 14 . 
The numerical data used here is the same as was used by 
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those researchers. These values are \* =1.5613, 0 f =i, and 
d 1 =d 2 =4, where d x and d 2 are the normalized dimensions of 
the quarter square section. 

Because the density is constant throughout the 
solution domain, no effects of flow in the liquid were 
considered and only the energy equation was solved. For 
this problem, a total of 49 equal square elements were 
used with a time step corresponding to a t of 0.005. 

Figure 5.4 shows the calculated interface at various times 
during the solidification. The interface locations are 
also presented in terms of fraction of solidified matter 
along the diagonal and at the insulated boundaries in 
Figures 5.5 and 5.6, respectively. The calculated results 
compare quite well with the results from other 
researchers. Though no exact solution to this problem 
exists, a similarity solution for an infinite medium is 
also given. The calculated results of the present 
analysis for the finite medium compare favorably with the 
infinite medium analysis initially. At later times, the 
infinite medium solution predicts a slightly faster 
solidification. This is expected since end effects in the 
infinite medium allow for heat conduction out of the 
corner , while ends for the finite medium are effectively 
insulated. 
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Figure 5.5 Calculated Interface location Along the 
Diagonal for Case 2 



Figure 5.6 Calculated Interface Location Along the 
Center Line for Case 2 
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Case 3: Buoyancy-Driven Convection 

Flow inside a square cavity is one of the simplest 
problems in convection and is often used to test the 
validity of fluid analysis methods. Case 3 is an example 
of buoyancy driven flow in a square cavity. Figure 5.7 
describes the solution domain and boundary conditions. 

For this case, one side of the square cavity is maintained 
at a constant temperature of T : . The opposite side is at 
a higher temperature T 2 . The velocities of the fluid are 
prescribed to zero at the container wall. The Prandtl 
number of the fluid was chosen to be 1. 

The dimensionless Rayleigh number, defined by 


Ra 


gp(T 2 -TJD 3 
va 


is used to characterize the flow. The coefficient of 
thermal expansion, 0, in the definition of the Rayleigh 
number above is defined by 


0 = 



1 

AT 


where the 0 subscript denotes some reference state, and AT 
is the change in temperature from that reference state. 
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Figure 5.7 Description of Case 3 


For this problem, a total of 49 equal sized square 
elements were used. Figures 5.8, 5.9, and 5.10 show the 
calculated fluid velocity fields and temperatures at 
Rayleigh numbers of 10 3 , 10* , and 10 5 . No cases were 
calculated for higher Rayleigh numbers, however, no 
stability problems were observed at a Rayleigh number of 
10 5 . The calculated results compare well qualitatively 
with similar analyses by Pepper and Cooper 49 . Pepper and 
Cooper also compiled data from the literature for the 
buoyancy driven cavity flow problem described above. This 
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Figure 5.8 Calculated Fluid Velocity Vectors and 
Normalized Temperature Contours for 
Case 3 at a Rayleigh Number of 10^ 



Figure 5.9 Calculated Fluid Velocity Vectors and 
Normalized Temperture Contours for 
Case 3 at a Rayleigh Number of 10* 
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data is presented in terms of an average Nusselt number 
given by 


Nu 



9T 

ax 


dy 


The Nusselt number above is dimensionless by normalizing 
with respect to the wall temperatures and letting the 
dimension D be 1. Results from the present analysis are 
given in Table 5.2 and graphically in Figure 5.11 along 
with the results from other researchers. 


Table 5.2 Calculated Nusselt Numbers for Case 3 


Ra 

Nu 

10 3 

1.11 

10 4 

2.40 

10 5 

5.17 


Figure 5.11 shows excellent agreement between the 
present analysis and previously published data. 

In addition to this problem, several other cases 
involving free and forced convection were studied during 
this research. One case considered driven cavity flow 
(forced convection) without buoyancy. In this case three 



Figure 5.10 Calculated Fluid Velocity Vectors and 
Normalized Temperature Contours for 
Case 3 at a Rayleigh Number of 10 5 



Rayleigh Number, Ra 

Figure 5.11 Average Nusselt Number Versus 
Rayleigh Number for Case 3 
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sides of the square cavity are maintained at a constant 

temperature of T 1 . The fourth side is at a temperature 

T 2 , and moves with a velocity V. The predicted velocity 

and temperature contours compared well qualitatively with 

13 

similar analyses by Chen, et al. 

Case 4: Combined Con duction, — Convection ,. 

and Phase Change 

The last case gives the solution to the general problem 
combining conduction, convection, and phase change. 

Figure 5.12 shows the solution domain along with the 
boundary and initial conditions. The problem consists of 
a material inside a container of square cross-section. 

The top and bottom of the container are maintained at 
constant temperatures above and below the fusion 
temperature of the material. Since the container is 
initially at rest and not subject to a gravitational 
field, there is no motion in the liquid phase. The 
initial temperature distribution varies linearly between 
T and T and flat phase front has formed perpendicular to 
the direction of the temperature gradient. At time zero, 
the material is subjected to an acceleration causing 
buoyancy-driven convection in the fluid. This convective 
flow changes the heat transfer and thus affects the 

location of the phase front. 

Table 5 . 3 gives the numerical values for the boundary 
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Figure 5.12 Description of Case 4 


conditions, material properties, and geometry of Case 4. 

The complexity of this problem is apparent from the 
number of physical constants listed in Table 5.3. Even if 
dimensional analysis were applied to this problem, over 
ten dimensionless groups would result. That many 
dimensionless numbers would not significantly improve the 
description of the problem, thus only the conventional 
dimensionless parameters are given here. 
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Table 5.3 Physical Properties and Conditions for Case 4 


Symbol 

Description 

Value 

T i 

Fusion Temperature 

0.0°C 

T i 

Prescribed Low Temperature 

-0 . 5 °C 

T 2 

Prescribed High Temperature 

1 . 5 °C 

p. 

Solid density 

1000 kg/m 3 

pj. 

Liquid density § T f 

1000 kg/m 3 

p 

Liquid thermal expansivity 

0.001 /°C 

p 

Dynamic viscosity 

0.001 kg/ 

c 

Solid specific heat 

1000 J/kg°C 

8 

C x 

Liquid specific heat 

1000 J/kg°C 

L 

Latent heat of fusion 

1000 J/kg 

k 

Solid thermal conductivity 

1 W/m°C 

K 

Liquid thermal conductivity 

1 W/m°C 

g 

Acceleration 

0.02 m/s 2 

D 

Length of container side 

0.1m 

t 

time 

s 


The Rayleigh number, characterizing the buoyancy 
driven convection, and given by 


Ra 


g/J(T 2 - T f )D 3 
va 


3x10* 


The Prandtl number, characterizing the fluid / s ratio 
of momentum and thermal diffusivity, and given by 


Pr = — = 1.0 
a 


The Jakob number, characterizing the material's ratio 
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of specific heat to latent heat capacity, and given by 


C (T - T ) 

Ja = — ! i 2— = o.5 


In addition, we may define dimensionless temperatures 
given by 



and 


0 

2 


T - T 

2 f 

T - T 

2 1 


0.75 


Finally, we may define a dimensionless time parameter 
given by 


a t 



For this problem a total of 49 equal square elements 
were used and the time step corresponded to a t of 0.0002. 
Figures 5.13 through 5.23 show the calculated results for 
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Case 4 at various times during the phase change process. 
These results include normalized isotherms, velocity 
vectors, and the location of the liquid/solid phase front. 
The development of the fluid flow and its effect on the 
heat transfer can be seen at the early time steps. The 
interesting aspect of this example is the definite 
influence of the fluid flow on the heat transfer and the 
resulting movement of the phase front. Steady state is 
reached at t =0.2 with a phase front significantly 
different from the flat front formed initially under 
conditions of no convective flow. 
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Figure 5.13 Normalized Isotherms, Velocity Vectors, 
and Phase Distribution for Case 4 
at t=0.0 







Figure 5.14 Normalized Isotherms, Velocity Vectors, 
and Phase Distribution for Case 4 
at t=0. 001 
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Figure 5.15 Normalized Isotherms, Velocity Vectors, 
and Phase Distribution for Case 4 
at t=0 . 005 
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Velocity Vectors, 
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Figure 5.17 Normalized Isotherms, Velocity Vectors, 
and Phase Distribution for Case 4 
at t=0 . 02 





Figure 5.18 Normalized Isotherms, Velocity Vectors, 
and Phase Distribution for Case 4 
at t=0 . 03 
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Nomalized Isotherms, Velocity Vectors, 
and Phase Distribution for Case 4 
at t=0 . 04 


Figure 5.19 
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Figure 5.20 Normalized Isotherms, Velocity Vectors, 
and Phase Distribution for Case 4 
at t=0.05 
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Figure 5.21 Normalized Isotherms, Velocity Vectors, 
and Phase Distribution for Case 4 
at t=0.1 
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Figure 5.22 Normalized Isotherms, Velocity Vectors, 
and Phase Distribution for Case 4 
at t=0. 15 


102 



Figure 5.23 Normalized Isotherms, Velocity Vectors, 
and Phase Distribution for Case 4 
at t—0.2 




103 


Computer Resource Usage 

Computer usage is of importance to most numerical 
modelers. Table 5.4 shows the computer time usage for 
each of the cases presented. The computer system used was 
an AMDAHL 5870 with an MVSXA operating system. It is 
important to note, that these cases were not fully 
optimized in terms of domain discretization or time step 
to provide minimum CPU times. 

Table 5.4 CPU Times of Verification Cases 


Case 

CPU Time, 
seconds 

CPU Time 

# Time Steps 

1 

696 

1.74 

2 

16470 

8.24 

3 

15035 

150.4 

4 

145700 

145.7 


The CPU times given in Table 5.4 are long and could 
result in significant expense on a pay for time computer 
system. The intended use of this analytical model, 
however, is in the aerospace industry and government with 
institutional computational facilities devoted to such 
tasks. Where no other similar analysis tool is available. 
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the comparison may be between these computer costs and the 
need for a space flight experiment, for example, the cost 
of which can also be substantial. 

It should also be noted that these CPU times are not 
unusual for modeling phase change problems with 
convection. Schneider 52 cites researchers quoting CPU 
times of 50000 seconds on a CDC 6500 computer, for similar 
phase change problems. Further discussion on computer 
usage and possible areas of improvement can be found in 
the next chapter. 



CHAPTER 6 


CONCLUSIONS AND RECOMMENDATIONS 


The work presented develops a solution approach to 
the co mb ined conservation laws of energy, momentum, and 
mass. This approach is quite general and provides a 
method to analyze a variety of fluid and thermal problems 

including those with phase change. 

The finite element method, being an integral method 
provides a natural means to implementing the internal 
energy formulation of the phase change problem. Using 
elements with quadratic interpolation functions, the 
interface can be tracked quite accurately. Finite 
difference formulations to date only provide information 
on the interface location to within on mesh spacing. 
Results from those analyses show jagged interfaces (see 
Schneider 52 ) that can influence convective flow in the 
liquid. 

The analytic approach was implemented in a computer 
program and results were verified by investigating 
individual phenomenon and comparing with known solutions. 
In general, the approach yields solutions with good 
engineering accuracy. The predicted results for a problem 
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similar to that described in Chapter 3 were presented for 
comparison by other researchers who may develop similar 
analytical methods in the future. 

One area of concern for the present approach is the 
high computer resource usage. The majority of the 
computer time is incurred in forming and then inverting 
the large set of system equations needed to converge the 
nonlinearities and achieve a valid solution. 

One way of reducing the computer time requirements by 
improving the algorithm may be to incorporate a Newton- 
Raphson iteration scheme to converge for the 
nonlinearities. For a further discussion on the 
application of of such methods, see Geradin, et al. 25 . 

Another approach to reducing computational times 
would be to improve the simultaneous equation solver. The 
matrices representing the system equations are 
characteristically sparse with nonzero coefficients 
located close to the diagonal. The use of a banded matrix 
solver could significantly reduce computer time. 

Fortunately, future trends will continue to reduce 
the computational times and costs. The operation of the 
computer program on a 3090 class computer with vector 
optimizing hardware should reduce CPU times by about an 
order of magnitude. This trend of increased computational 
capability of the hardware will continue in the future. 
Also the development of efficient large matrix solvers is 
undergoing much research and is a very important aspect of 
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numerical modeling of fluid/thermal systems. 

There are several obvious extensions to this work. 
First, this approach can easily be modified to analyze 3- 
dimensional space systems. Nothing in the formulation 
should prevent this and only the available computational 
capabilities might present a restriction on the geometric 
complexity of the problem. 

Further investigation is also warranted in the use of 
consistent mass matrices for the inertia terms in the 
energy and momentum equations. Though the lumped mass 
matrices are by far the most common, there is concern that 
the method of lumping might contribute to inaccuracies 
particularly for highly distorted curve-sided elements. 

Though I personally believe that the internal energy 
or enthalpy method is the only practical method for 
modeling the phase change problem in 3 -dimensional space, 
development of other methods is warranted. The work of 
Chang and Brown ^ is particularly interesting, although 
their moving boundary -moving mesh techniques would be very 
difficult to implement where the thermodynamic phases are 
fragmented. 

For application to the space environment, other 
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forces not addressed here can become important. Siegel 
provides a good overview on the effects of reduced gravity 
on heat transfer. The incorporation of a surface free 
energy model would be especially useful, particularly for 
very low gravity conditions. Pearson^®, Labus^^, and 
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Merte 40 provide interesting approaches to handling the 
free surface problem, however, application to the finite 
element method remains fertile ground for research. 

Finally, the incorporation of a turbulence model into 
the analysis could significantly extend its application. 
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APPENDIX A 


COMPUTER PROGRAMMING APPROACH AND DOCUMENTATION 


This appendix discusses the general programming 
approach and implementation used in this research. 

Included is a description of the program functions and 
important global variables. This is followed by a 
flowchart of the major functions. 

The programming language APL was chosen for its 
reduced programming development time, its inherent matrix 
manipulation capabilities, and because it is easily 
transported to various computer hardware systems. APL was 
originally developed by Iverson 33 as a- general 
mathematical notation and later implemented as a computer 
programming language. A detailed description of the 
programming language APL is given by Gilman and Rose 26 . 

APL has several significant advantages over more 
"conventional" languages such as FORTRAN, PASCAL, etc. 
Because APL is a symbolic vector language the source code 
is typically at least two to three times shorter than most 
other languages and can be developed about four to ten 
times faster. Table A.l shows a comparison between 
FORTRAN and APL programs to produce the sum of all the 
numbers greater than 50 in a set of real numbers. The APL 
program is significantly shorter. Several characteristics 
contribute to APL's brevity. APL processes aggregates of 
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data, or arrays, with much the same syntax as single 
numbers, and therefore essentially eliminates the need for 
structures such as loops. Data is maintained unformatted 
within the APL environment, so that no logical ties to the 
host operating system files is required. In addition, all 
numbers are stored and operations performed using double 
precision. 


Table A. 1 FORTRAN and APL Programs to Sum the Real 

Numbers In a Set That are Greater Than 50 


FORTRAN 

APL 

DOUBLE PRECISION SUM,X 
READ (5, 100) LIMIT 
100 FORMAT (110) 

SUM=0 . 0D0 

200 DO 400 1=1, LIMIT 
READ(5, 300) X 
300 FORMAT (F10. 5) 

IF (X. LE. 50 . ) GO TO 200 
400 SUM=SUM+X 

WRITE ( 6 , 500) SUM 
500 FORMAT (F20. 5) 

STOP 

END 

+/(X>50)/X 


Most APL systems today use an interpreter though some 
compilers are available. The interpreted versions lend 
themselves very well to program development because the 
environment is interactive, no compiling is required, and 
the debugging facilities work directly with the source 
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code. Compiled versions however can run significantly 
faster for some problems (particularly those that are 
highly iterative) . For the problem presented in this 
research, however, the computer time usage is dominated by 
inversion of the large system matrix. This inversion is 
performed by a highly efficient APL primitive function and 
would not benefit significantly from using a compiled 
external function. This APL matrix inversion function 
will also take advantage of super computer class vector 
processing hardware, when available, to further improve 
computational performance. 

The following pages contain descriptions of the 
program functions and important global variables. This is 
followed by the calling structure (flow chart) of the 
major functions. 
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FUNCTIONS 


ADJS adjusts time step based on ability to converge and field values 

INTERATIVE SOLUTION ( FOR NON-LINEARITIES ) AT EACH TIME STEP 
AXNTCON 

INTEGRATION CONSTANTS FOR USE WITH AINTGRT 

AREA 

CALCULATES AREAS AND VOLUMES OF ELEMENTS 

BDY 

EVALUATES R VECTOR FROM BOUNDARY AND INTERNAL CONDITIONS 
BDY2 MODIFIES ELEMENT R MATRIX FOR BOUNDARY COND. TYPE 2 ( PRESCRIBED FLUX ) 
BDY2 SUB 8 - FUNCTION OF BDY2 TO PRESCRIBE EACH ( ELEMENT , SIDE COMBINATION) 

CHKC CHECKS CONVERGENCE OF EF , UF , VF AND PF WITH LAST ITERATION VALUES 

CHKCNVLS 1 , „ 

CHECKS FIELD VALUES FOR CONVERGENCE 

CHKEF CHECKS ENERCY FIELD VALUES TO SEE IP WITHIN PROPERTY DATA RANCE 

CBKI CHECK PRESRCRIBED PRESSURE NODE " PRNODE " SPECIFIED .IN INPUT 
CHKINPUT 

CHECK INPUT PARAMETERS 
CHKSS 

CHECK IP STEADY STATE HAS BEEN REACHED 
CHKTF 

CHECKS TEMPERATURE FIELD VALUES IF WITHIN PRESCRIBED VALUES 

CLEA ERASES VARIABLES NAMED IN CLV EXCEPT THOSE WITH ALT. CHARACTER NAMES 
COORXE 

CENERATES XI-ETA COORDINATES OF THE NODES 
CPUCHK 

RETURNS A 0 IF CPU TIME LIMIT IS EXCEEDED 
CPUTIME 

RETURNS CPU SECONDS USED SINCE TSTART - WAS ISSUED 

DFN CALCULATES DERIVATIVES OP FIELD VARIABLE AT THE ELEMENT NODES 
DIAG 

FORMS DIAGONAL MATRIX FROM A VECTOR X 
EGYASN 

ASSIGN ENERGY FIELD VALUES 
EGYBAR 

FORMS ENERGY EQUATIONS FOR SOLVING TRANSIENT RESPONSE 

ECYEVL 

EVALUATES ENERGY MATRICES 
E GY MAT 

CONSTRUCT THE ENERCY ELEMENT MATRICES 
ENERGY 

FORMS AND SOLVES THE TRANSIENT ENERGY EQUATION 
FLDMATLCM 

FORMS CM MATRIX FOR FLDMAT 
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FLDMAT&KA 

FORMS KA MATRIX FOR FLDMAT 
FLDMAT LKAE 

FORMS KAE MATRIX FOR FLDMAT 
FLDMAT LKNN 

FORMS K 11 K 12 K21 K 22 MATRICIES FOR FLDMAT 
FLDMAT LLCN 

FORM LC 1 AND LC 2 MATtf/CITO TO/? TOflMAr 
FLDMAT LLN 

FORMS LI AND L2 MATRICIES FOR FLDMAT 
FLDMATLMM 

FORMS MASS MATRIX FOR FLDMAT 
FLDMATLR 

FORM R VECTORS FOR FLDMAT 
FLDMATLRC 

FORMS RC VECTOR FOR FLDMAT 
FLDMAT LRP 

FORMS AND ASSIGNS RP MATRICIES FOR FLDMAT 
FLDSLV 

SOLVES EQUATIONS AND ASSIGNS FIELD VALUES 
FLDSLV LSTAR 

RETURNS FIELDS FROM LAST ITERATION 

FLOW 

FORM AND SOLVE THE FLUID FLOW EQUATIONS 
FLWASN 

ASSIGN FLUID FLOW FIELD VALUES 
FLWBAR 

FORMS FLOW EQUATIONS FOR SOLVING TRANSIENT RESPONSE 
FLWEVL 

EVALUATES FLUID FLOW ELEMENT MATRICES 
FLWMAT 

CONSTRUCT THE FLUID FLOW MATRICES 
FMFEQL 1 

FORM ENERGY FIELD EQUATIONS 
FMFEQL 2 

FORM FIELD EQUATIONS 
FRONTEND 

PERFORMS UPFRONT ONCE ONLY FUNCTIONS 
CRAW EC 

RETURNS GRAVITATIONAL BODY FORCE VECTOR 
GRIDLELN 

GENERATES QUADRATIC NODE NUMBERS ( ELEMENT BASIS ) FOR CRIDGEN 
GRIDLLNODE 

GENERATES LINEAR NODE NUMBERS ( ELEMENT BASIS ) FOR GRIDGEN 
GRIDLRSE 

GENERATES REGION SIDE ELEMENT NUMBERS 
CRIDLRSLN 

GENERATES REGION SIDE LINEAR NODE NUMBERS FOR GRIDGEN 
GRIDLRSN 

GENERATES REGION SIDE QUADRATIC NODE NUMBERS FOR GRIDGEN 
GRIDGEN 

GENERATES ELEMENT NODE NUMBERING 
INBDY 

SPECIFY REGION BOUNDARY CONDITIONS 
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INGEOM 

SPECIFY GEOMETRY INPUT PARAMETERS 
INITIAL 

INITIALIZE VARIABLES 
INITMP 

INITIALIZE MATERIAL PROPERTIES : CASE WHERE PRESSURE NODE PRESCRIBED 
INITMPROP 

INITIALIZE MATERIAL PROPERTIES USING AVERAGE CONDITIONS 
INITPRES 

INITIALIZE PRESSURE FIELD 
INITTEMP 

INITIALIZE TEMPERATURE 
INITVEL 

INITIALIZE FLUID VELOCITIES 
INPLYGN 

FINDS IF POINTS XY ARE IN POLYGONS Q 
INPROG 

SPECIFY PROGRAM OPERATION INPUT PARAMETERS 
INPUT 

SPECIFY AND PRINT OUT INPUT PARAMETERS 
IN REG 

SPECIFY INITIAL R AND Z COORDINATES OF REGION 
INTQR 

INTERPOLATES BY FINDING CLOSEST REGION AND USING DBL . QUAD. REGRESS . 
INTQR&REG 

CALLED BY INTQR , IT PERFORMS DOUBLE QUADRATIC REGRESSION BY REGION 
INTQRLUNN 

UNNORMALIZE A MATRIX "A" ( NESTED ARRAY ) WRT . RANGE 
JACCHK 

CHECK IF DETERMINANT OF JACOBIAN HAS A SIGN REVERSAL 
LSHAPE 

CALCULATES LINEAR SHAPE FUNCTIONS AND THEIR DERIVATIVES 

LUMP 

LUMPS THE CAPACITANCE MATRIX BY ROWWISE SUMMATION 
LX 

LATENT EXPRESSION FOR PHASTRAN 

MAP 

MAP XI-ETA COORDINATES OF NODES INTO X-Y SYSTEM 
PHASTRAN 

MAIN CONTROL FUNCTION FOR PHASE CHANGE ANALYSIS MODEL 
PRSCRB 

PRESCRIBES FINITE ELEMENT EQUATIONS IN THE GLOBAL VARS. KBAR AND RBAR 
PRSIDE 

MODIFIES PFEQ WHICH DEFINES PRESCRIBED BOUNDARY CONDTIONS 
PRSIDELADJ 

ADJUSTS ORIGINAL NODE POSITIONS IN A MATRIX FOR FIELD TYPE 
PRSIDELFLD 

PRESCRIBES BOUNDARY CONDITIONS FOR A FIELD 
PRSIDE LTEMP 

CONVERTS TEMPERATURE BOUNDARY CONDITIONS TO ENERGY BC % S 
PRSPRES 

PRESCRIBES PRESSURE NODEM 
PRSP 1 

PRESCRIBES PRESSURE NODEM: CASE WHERE TOTAL MASS CONSTRAINED 
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PRSP 2 

PRESCRIBES PRESSURE NO DEM: CASE OP OPEN SI STEM 
PRSVLC 

RETURNS LOCATIONS OF PRESCRIBED VALUES FOR THESE FIELD TYPES 
QLGSHP 

CALCS. QUADRATIC LACRANCIAN SHAPE FUNCTIONS AND DERIVATIVES 
QUADRCXY 

PERFORMS QUADRATIC REGRESSION ANALYSIS IN 3 DIMENSIONS X Y Z 
RCOND 

RETURNS SUB -VECTOR FOR INTERNAL CONDUCTION 
RDUPL 

RETURNS A BOOLEAN FOR REDUCING DUPLICATE VALUES LEAVING ONLY THE 1ST 
REDUCE 

REDUCES (BY SUMMATION) A VECTOR WITH MULTIPLE INDICIES 
REMDUPEL 

REMOVES DUPLICATE ELEMENTS OF X 
RESULTS 

DISPLAYS FIELD VARIABLES 
RHOZERO 

CALCULATES THE REFERENCE DENSITY 
RPRES 

FORMS PRESSURE RESULTANT VECTOR 
SAVEFIELDS 

SAVES FIELDS AT EVERY DFN TIME STEP , PUTS IN NESTED ARRAY SAVEF£,D$ 
SETLAST 

SETS VARIABLES FROM LAST TIME ITERATION 
SETSOLID 

SETS VELOCITIES TO ZERO FOR NODES IN SOLID STATE 
SETSTAR 

SETS EF> UF , VF AND PF AT LAST CONVERGENCE ITERATION 
SIFAC 

CALCULATES THE SIDE INTEGRATION FACTOR FOR USE WITH SINTGRT 
SINTCON 

INTEGRATION CONSTANTS FOR USE WITH SINTGRT 
STATPRES 

CALCULATES THE STATIC PRESSURE DISTRIBUTION 
STRIPE 

CREATES STRIPE LINE BORDER 
TIMECHK 

RETURNS A 0 IP TIME LIMIT IS EXCEEDED 
TIMESTEP 

TIME STEPPING FUNCTION 
UPPROP 

UPDATES PROPERTIES ON GLOBAL NODE BASIS 
UPPROP LTEMP 

PLACES PRESCRIBED TEMPERATURES IN UPDATED PROPERTIES 
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VARIABLES 


ANI ORDER OR CAUSS-LECENDRE QUADRATURE FOR AREA INTEGRATION 

AXEI X I AND ETA COORD . USED IN AREA NUM. INTEGRATION 
BC 

MODIFIED BOUNDARY CONDITIONS OF REGIONS 

CEQ 

CAPACITANCE MATRIX FOR FLOW EQUATIONS 
QI 

~ CONSTANTS FOR FUNCTION INTGRT 
CIL 

CONVERGENCE ITERATION LIMIT 

CM AT 

HEAT CAPACITANCE MATRIX 
CPULIM 

LIMIT FOR CPU TIME 

CSL 

CONSTANTS FOR LINEAR SHAPE FUNCTIONS 
ONSTANTS FOR QUADRATIC LAGRANGIAN SHAPE FUNCTIONS 

dxt 

DER. SHAPE FNCS, WRT X FOR ELMNT . TYPE (NESTED ARRAY ) 

DYT 

DER. SHAPE FNCS. WRT Y FOR ELMNT. TYPE (.NESTED ARRAY ) 

EC 

ENERGY CALCULATION CONTROL (0 -NO THERMAL CALCS.) 

EF 

INTERNAL ENERGY 
EFINIT 

INITIAL INTERNAL ENERGY 
EFSTAR 

INTERNAL ENERGY OF LAST CONVERGENCE ITERATION 

ELN 

NODE NUMBERS FOR EACH QUADRATIC ELEMENT 

EREQ 

RESULTANT VECTOR FOR ENERGY EQUATIONS 

ERR 

ALLOWABLE FIELD CONVERGENCE ERROR 

FELT 

FIELD ELEMENT TYPES ( 1 -LINEAR , 2-QUAD) 

FT 

FIELD TYPE (1-INTERNAL ENERGY , ETC ) 

FC 

FLUID CALCULATION CONTROL (0 -NO FLOW, 1-FLOW CALCULATED) 

FREQ 

RESULTANT VECTOR FOR FLOW EQUATIONS 
GF 

GEOMETRIC FACTOR FOR X-Y OR R-Z COORDINATE SYSTEMS 

GRZ 

GRAVITATIONAL CONSTANTS IN X AND Y DIRECTIONS 

ICTL 

ITERATION CONTROL (0 -SUCCESSIVE SUB. 1 -NEWTON-RAPHSON) 
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KBAR 

MODIFIED STIFFNESS MATRIX FOR TRANSIENT FLOW EQUATIONS 

KEQ 

STIFFNESS MATRIX FOR FLOW EQUATIONS 
KT 

THERMAL CONDUCTIVITY 

LNODE 

NODE NUMBERS FOR EACH LINEAR ELEMENT 

LSS 

SHAPE FUNCTIONS FOR SIDE INTEGRATION 
ND 

NUMBER OF DIVISIONS PER SIDE PER REGION 
NE 

NUMBER OF ELEMENTS PER REGION 

NST 

SHAPE FUNCTIONS FOR ELEMENT TYPES (.NESTED ARRAY) 

PF 

INTERPOLATED PRESSURE FIELD OF FLUID 

PFEQ 

MATRIX POS. AND PRSCRB . VALUES OF FLUID 
PFINIT 

INITIAL PRESSURE FIELD 
PFSTAR 

PRESSURE OF LAST CONVERGENCE ITERATION 

EPHI 

ENERGY FIELD VARIABLE SOLUTION 

FPHI 

FLOW FIELD VARIABLE SOLUTION 
PRNODE 

PRESCRIBED PRESSURE NODE INFORMATION 

QAIF 

QUAD. AREA INTEGRATION FACTOR 

QSIF 

SIDE INTEGRATION FACTOR 

QSS 

SHAPE FUNCTIONS FOR SIDE INTEGRATION 

RBAR 

MODIFIED RESULTANT VECTOR FOR TRANSIENT FLOW EQUATIONS 

REQ 

RESULTANT VECTOR FOR FLOW EQUATIONS 

RHO 

DENSITY OF MATERIAL 
RE LAST 

LAST TIME ITERATION RESULTANT VECTOR FOR ENERGY 
RFLAST 

LAST TIME ITERATION RESULTANT VECTOR FOR FLOW 
RN 

REGION NUMBER 
RSELMT 

REGION SIDE ELEMENTS 

RSLN 

REGION SIDE LINEAR NODES 

RSN 


REGION SIDE NODES 
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RZC 

X AND Y COORDINATES OF CONTAINER 

RZE 

X AND Y COORDINATES OF ELEMENTS 

RZN 

X AND Y COORDINATES OF NODES 

RZR 

X AND Y COORDINATES OP NODES OF REGIONS 
SAREA 

SIDE AREAS OP RING ELEMENTS 

SSC 

STEADY STATE CONTROL (1 -EXIT EARLY IF REACH STEADY STATE ) 
SF 

STATE (OF THE MATERIAL) FIELD 

SFN 

SAVE FIELDS EVERY SFN TIME STEPS 

SNI 

ORDER OR GAUSS-LEGENDRE QUADRATURE FOR SIDE INTEGRATION 
STATMSG 

STATUS MESSAGES 

SXEI 

XI AND ETA COORD. USED IN SIDE NUM. INTEGRATION 

TCTL 

TIME START , END AND INCREMENT CONTROL 
TF 

TEMPERATURE FIELD 
THETA 

TRANSIENT ALGORITHM CONTROL PARAMETER 

TIME 

RECORD OF TIMES FGR EACH TIME STEP 
TIMELIM 

TIME LIMIT ON RUN 
UNIT 

INITIAL TEMPERATURE 

TSC 

TIME STEP CONTROL ( 0 -CONSTANT , 1 -VARIABLE ) 

UP 

U VELOCITY FIELD 
UFINIT 

INITIAL U VELOCITY FIELD 
UFSTAR 

U VELOCITY FIELD OF LAST CONVERGENCE ITERATION 
VF 

V VELOCITY FIELD 
VFINIT 

INITIAL V VELOCITY FIELD 
VFSTAR 

V VELOCITY FIELD OF LAST CONVERGENCE ITERATION 

VIS 

FLUID VISCOSITY 

VOL 

VOLUME OF THE RING ELEMENTS 

XEN 

XI AND ETA COORD, OF NODES URT REGION COORDINATES 


FLOWCHART OF WORKSPACE 


LX 

:CLS 

: STRIPE 
PH ASTRAS 

: FROSTEND 

: CLEANUP 
: INPUT 

UNCEOM 

: 1NREC 
UNBDY 
: INPROC 
: CRIDGEN 

: GRIDLRSN 
: CRIDLRSE 
: GRIDLRSLN 
: GRIDLELN 
: GRIDLLNODE 
: COORXE 

■.MAP 

■.INITIAL 

: CHKINPUT 

: CHKINPRES 

•.A I NT CON 

: LSHAPE 
: QLGSHP 
: JACCHK 
: SINTCON 

: QLCSHP 

-.LSHAPE 

: SI FAC 

-.AREA 
: INITTEMP 

: INITVEL 
: INITPRES 

lINITMP 

: INTQR 

: INPLTGN 

: INTQRLREG 

: REMDUPEL 

-.QUADRCXY 

: INTQRLUNN 
: UPPROP 

: INTQR 
'.INPLYCN 

: INTQRLREC 

: REMDUPEL 
: QUADRCXY 
: INTQRLUNN 
: UPPROPLTEMP 

: ST ATP RES 

: INITMPROP 
: RPRES 
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:DFN 

: QLGSHP 

: DFN 

: QLGSHP 

zRCOND 
: GRAVVEC 

iRHOZERO 

: INTQR 

: INPLYGN 
: INTQRLREG 

zREMDUPEL 
: QUADRGXY 
: INTQRLUNN 

: SETLAST 
: SETSTAR 

TIMESTEP 

: CP£/Ctf K 

: CPUTIME 
: r IMPCPK 
:ACAIJV 

: SETSTAR 
: ENERGY 

: EGYEVL 

lEGYMAT 

iFLDMATLCM 
zFLDMATLRC 
: BDY 

\BDY2 

: BDY2LSUB 

: RCOND 
: FLDMATLKAE 
: FLDMATLRP 
: RPRES 

zDFN 

: (JLCSPF 

IDFN 

: QLGSHP 

iFMFEQL 1 
jFCYSAF 

:Li/MP 

iPIAG 

jPFSIDF 

:P/?SI0FAF£D 

rPPSJDPATPMP 
: INTQR 

: IJVPLrCN 
: INTQRLREG 

iREMDUPEL 
: QUADRGXY 

IPRSIDELADJ 
iRDUPL 
: PRSCRB 

iRDUPL 
: PRSVLC 
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zFLDSLV 

: FLDSLV&STAR 
: EGYASN 

z FLOW 

: FLWEVL 

ZFLWMAT 

z FLDMATLMM 
zFLDHAT&KA 
zFLDMAT&KNN 
: FLDMATLLN 
: FLDMATLLCN 
zFLDMAT&R 

: CRAW EC 

z RROZERO 

zINTQR 

: INPLYGN 
: INTQRbREG 

z REMDUPEL 
IQUADRGXY 
: INTQR&UNN 

: BDY 

zBDY 2 

ZBDY2LSUB 

: BDY 


ZBDY 2 

zBDYibSUB 

z FMFEQ&2 
zFLWBAR 

zLUMP 

zDIAG 

zPRSIDE 

zPRSIDE&FLD 

z PRSIDEbTEMP 
z 

: INPLYGN 
zINTQRLREG 

z REMDUPEL 
: QUADRGXY 
z INTQRL UNN 
zPRSIDE LADJ 
zRDUPL 
zPRSPRES 

z PRSP1 
z PRSP2 
iPRSCRB 

zRDUPL 

zPRSVLC 

zFLDSLV 

zFLDSLV LSTAR 
z FLWASN 
z UPPROP 

zINTQR 

z INPLYGN 
zINTQR tREG 

Z REMDUPEL 
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IQUADRGXY 
-.INTQRbUNN 
: UPPROPLTEMP 
: CRKCNV 

: CRKCNV US 1 
: CRKCNV LS 1 
sCPKCAWASl 
:CPKC/WAS1 
:SEr50Lii? 
sprso^ip 
ADJSTEP 

: CffKFF 
: CRETE 
: INITIAL 

: CtfKIJVPC/r 

: CHKINPRES 
lAINTCON 

:LSHAPE 
: QLGSHP 
:JACCRK 

: s r/vrow 

:£StfdP£ 
t SI FAC 

lAREA 
: INITTEMP 
: INITVEL 
: INITPRES 

: IN ITMP 
: 

:IJVPZ,ycW 

: 

IFEWWPPX 
: GiMMCXy 
iintqrlunn 

: UPPROP 

: IJVrCP 

:i/vpz;y^ 

: INTQRLREG 

iREMDUPEL 
IQUADRGXY 
lINTQR&UNN 
: UPPROP&TEMP 

: STATPRES 
zINITMPROP 
:RPRES 

: DFN 

:QLGSRP 

: DFN 

: QLGSRP 

iRCOND 
: GRAVVEC 

iRHOZERO 

: INTQR 

lINPLYGN 
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xINTQRLREG . 

iREMDUPEL 

iQUADRCXI 


XINTQR&UNN 

iSETLAST 
: SET STAR 

: SETLAST 
-.SAVE FI ELDS 
l CHKSS 
: RESULTS 
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APPENDIX B 


COMPUTER PROGRAM LISTING AND RESULTS 


This appendix contains the APL source program listing 
used in this research. The input data for Case 4 
described in Chapter 5 is contained in the functions and 
global variables listed in this appendix. To run the 
program after loading the APL workspace enter the 
following: 

PHASTRAN # 

where # is the the number of divisions of a side of the 
solution domain. A value of seven results in 49 equal 
elements for this case. The results of such a run are 
given in the global variable RESULTS at the end of this 
appendix. 
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SYSTEM VARIABLE SETTINGS 
OCr: 1£~13 
DEC: . ,*0_" 

□ HT: 

□JO: 1 
QLC: 

□IX: 

DPP : 5 
DPR-. 

DPW : 79 
DRL: 16807 
DTZ: 0 

□J/4: 31226812 
0 NLT: 


ADDROWS 

[0] TRIABLE ADDROWS ROW-, l 

[1] a ADDS ROW(.S ) TO A TABLE FILLINC WITH BLANKS OR 0'S 


[2] R 

[3] TABLE*TOMATRIX TABLE a 

[4] ROW* TOMA TRI X ROW a 

[5] L*Cl+pTABLE)r~l + DROW a 

[6] TABLE*<(l+pTABLE) , DATABLE R 

[7] r«-MBIE,[l3((l + ptfOfcO,I)f/?OV a 


MAKE SURE TABLE IS A MATRIX 
CONVERT ROW TO A MATRIX 
FIND LARGEST OF COLUMNS 
RESHAPE- TABLE 
ADD ROW IS) 


ADJSTEP 

[03 CIC ADJSTEP NIT: MSG 

[1] R ADJUSTS TIME STEP BASED ON ABILITY TO CONVERGE AND FIELD VALUES 


[23 a 

[33 +<rsc=o)/o R 

[43 *(NIT>CIL)/QEC R 
[53 ♦ ( ~CHKEF ) /DEC a 

[63 *(~CHKTP)/[)EC R 

[73 -*-0 fl 

[83 ££C : TIME*TCTL [13 « 

[93 INITIAL R 

[103 NTS * 1 a 

[113 CIC * lOpO a 

[123 TCTLL3)*TCTLZ31*2 a 

[133 £ND:MSG*'TIME STEP CHANCED TO' a 

[143 STATUS MSG TCTLl 33 a 

[153 -*0 a 

[163 STOP: STATUS 'NOT CONVERGED ' a 
[173 ♦ A 


EXIT IF TIME STEP CAN'T CHANCE 
CHECK FOR TOO MANY ITERATIONS 
IF ENERGY OUT OF RANCE , DECREASE 
IP TEMP. OUT OF RANGE, DECREASE 
EXIT, TIME STEP IS OK 
SET BACK TIME TO BEGINNING 
REINITIALIZE PROBLEM 
DECREASE TIME STEP COUNTER 
REINITIALIZE ITER. COUNTER 
DIVIDE TIME STEP BY 2 
MESSAGE TO USER 
DISPLAY AND RECORD MESS ACE 
EXIT 

MESSAGE TO USER 
END EXECUTION 
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A I NT CON 
[ 0 ] 

[ 1 ] « 
[2] o 
[33 

[4] 

[5] 

[ 6 ] 

[7] 

[83 

[93 

[103 

[113 

[123 

[133 

[143 

[153 

[163 

[173 

[183 

[193 

[203 

[213 

[ 22 ] 

[23] 

[24] 


AINTCON N\B1\B2\D\N1\ETA\QJAC\NSHP\LAS\QAS 
INTEGRATION CONSTANTS FOR USE WITH AINTCRT 


Nl+l+ANI+N p 

ETA+(N1 f Nl)pNl+(QI[l; ; 3 ) LAN I ; 3 p 

AXEI+ ( 2 , {N1*N \ ) )p ( , *ETA ) , ( , ETA ) p 

LAS+LSHAPE AXEI p 

QAS+QLGSHP AXEI p 

QJAC+RZE JACOB QAS o 

QAIF+MDET QJAC P 

JACCHK QAIF P 

Bl4-P2P[; ; 1 3 5 73 BP LAS p 

B2+RZE BF QAS p 

tfSflP<-4 2 1 3*(l,Are,pAfSffP)pAtfffP+MS[l; 

p 

NSHP+ 4 2 1 3*(1 ,NE,pNSHP)pNSHP+QASll; 
NST+NST 9 cNSHP p 
04-4 1 2 3$(l,p£>)p£>4-£l[; ;1;3 p 
DXT^ cD p 

£>*■4 l 2 3$<1 •p0)p0*£lCl 52:3 p 
DYT+^d p 

£>4-4 1 2 3$( 1 » p0 )p04-E2 [ ; ; 1 ; 3 p 

£)Xr4-PXr i c£i p 

£> 4-4 l 2 3$ ( 1 , pD )pD+B2 [ ; ; 2 ; ] p 
pyr+pyr.cD « 


0P0BP OF GAUSS-LEGENDRE 
ETA COORDS. 

XI AND ETA COORDS . 

LINEAR SHAPE FUNCTION 
QUAD. LAGRANGIAN SHAPE FNS. 
JACOBIAN FOR QUAD. ELEMENTS 
QUAD . AREA INTEGRATE FACTOR 
CHECK JACOBIAN 
LINEAR GRADIANT MATRIX 
QUAD. GRADIANT MATRIX 
; 3 p SHAPE FNS. ELMNT. TYPE 1 
PUT IN GLOBAL VARIABLE 
• 3 p SHAPE FNS. ELMNT. TYPE 2 
ADD TO GLOBAL VARIABLE 
DERIV . WRT. X ELMNT . TYPE 1 
PUT IN GLOBAL VARIABLE 
DERIV. WRT. Y ELMNT. TYPE 1 
ADD TO GLOBAL VARIABLE 
DERIV. WRT. X ELMNT. TYPE 2 
ADD TO GLOBAL VARIABLE 
DERIV. WRT. Y ELMNT. TYPE 2 
ADD TO GLOBAL VARIABLE 
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AINTGRT 

CO] I+AINTGRT A;WE;WEX;N1 

[I] a INTEGRATES FUNCTION OVER AREA OF ELEMENT 
[23 a 

[33 Nl+ANI+1 a 

[4] +(lep ,A)/SCALAR a 

[5] +(*/(pQAIF)=2+qA)/&XT a 

[63 STATUS £3fiflSCCl3 a 

[7] +0 A 

[83 ffXr:4^x(2^ippi4)^(2<t>p4)pWJF a 

[9] +QALC a 

[10] SCALAR: A+AxQAIF a 

[II] QALC:WE+(N\ t N\)pN\*(Qll2\ ;] )ZANI;1 a 

[12] WEX+(pA)p*((*/(l+pA)),{l + f>A))Q(,WE)x 9 WE a 

[13] I++/lllWEX*A a 


IN XI -ETA COOR. SYSTEM 

1 +0RDER OF INTEGRATION 
CHECK IF A IS SCALAR 
CHECK IF p A IS LIKE p QAIF 
MESSAGE TO USER 
EXIT 

MULT . BY INTGRT. FACTOR 
JUMP TO FINISH INTEGRATION 

mult . fly lArrcflr. racroj? 

WEIGHTING FACTORS 
RESHAPE WEIGHT FACTORS 
NUMERICAL INTEGRATION 


AREA 

[0] AREAxRN 

[1] A CALCULATES AREAS AND VOLUMES OF ELEMENTS 

[2] A 

[3] SAREA+GFxSINTCRT 1 a SIDE AREAS OF ELEMENTS 

[4] VOL+AINTGRT 1 a VOLUME OF ELEMENTS 


BDY 

[03 

[13 

[ 2 ] 

[33 

[4] 

[5] 

[ 6 ] 
[7] 
[83 
[93 
[ 10 ] 
[ 11 ] 
C123 
[133 

[14] 

[15] 
[163 

[17] 

[18] 
[193 
[203 
[213 
[223 
[233 
[24] 
[253 
[26] 
[27] 


R+BDY FT ;B;M; NC ;N;D;S 
a EVALUATES R VECTOR FROM BOUNDARY 

A 

B+FT FSTCM BC a 
JVOBCsl] a 

R+( 4, (SNI+1 ) # iV£)pO a 

A 

LO0P:+(0=pNC)/NXTl a 
N+l+NC a 
NC+(N*NC)/NC a 
+( tf=l)/£0OP a 
B+N FSTCM B a 
D+RN FSTCM B A 
♦<O=x/p0)/£Xri a 
♦ (N- 2 )/g2 A 

M+ 1 BOUNDARY CONDITION NOT DEFINED 
STATUS M a 
+L00P a 

$2 :R+D BDY 2 R a 
+L00P a 

A 

NXTl:+(l 2-FELTIFT1 ) /LINEAR .QUAD a 
STOP a 

QUAD-.S+QSSL 1 5 5 3 a 
+NXT 2 a 

LINEAR: S+LSSZ 1;;3 a 
NXT2:R+* 2 3 l*(("l*pS),p*)p/? a 
R + 3 1 2 4$((p R)pS)*R a 


am? xwrflfliVAz; conditions 

BC'S FOR THIS FIELD TYPE 
LIST OF BC TYPES 
INITIALIZE R VECTOR TO ZERO 

LOOP ON BOUNDARY CONDITIONS 
TAKE FIRST BC TYPE 
REMOVE THIS TYPE FROM LIST 
IGNORE PRESCRIBED BC'S 
BC'S FOR THIS BC TYPE 
BC'S FOR THIS REGION NUMBER 
EXIT IF NO SUCH BC'S 
BOUNDARY CONDITION TYPE 2 
(BDY )' A WARNING TO USER 

DISPLAY AND RECORD MESSAGE 
CONTINUE LOOPING 
PRESCRIBED FLUX BC 
END OF LOOP 

CHECK IF LINEAR ELEMENTS 
WRONG ELEMENT TYPE 
QUADRATIC SHAPE FUNCTIONS 
FINISH CALC . OF R VECTOR 
LINEAR SHAPE FUNCTIONS 
RESHAPE R 

AND MULTIPLY BY SHAPE FNS . 
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BDY 2 

Cl] MODIFIES ELEMENT R MATRIX FOR BOUNDARY COND. TYPE 2 < PRESCRIBED FLUX) 
[2] p ; 1] IS THE SIDE NUMBERS 

C3] P DC; 2] IS THE PRESCRIBED FLUX VALUES 

C4] P 

C5] +(0=l*pD)/0 n 

C6] C+R o 
C7] Dl+ C C2]D p 
C8] BDY2tSUB”D\ p 
C9] C+R fl 


EXIT IF NO BOUNDARY CONDITIONS 
INITIALIZE RETURN VARIABLE 
NESTED ARRAY OF BOUNDARY COND. 
HANDLE EACH BOUNDARY CONDITION 
RETURN R VECTOR FROM BDYLSUB 


BDY2LSUB 


CO] 

Cl] 

C2] 

C3] 

C4] 

C5] 

C6] 

C7] 

C8] 


p SUB-FUNCTION OF BDY2 TO PRESCRIBE EACH (.ELEMENT , SIDE COMBINATION) 
Dt IS A NESTED ARRAY OF PRESCRIBED SIDE BOUNDARY CONDITIONS 


R IS THE R MATRIX FROM BDY 2 WHICH IS MODIFIED BY THIS FUNCTION 


S«-D1 Cl] p 

EL*-RSELMT(S‘, ) p 

PV*- ( (pR ) C2 ] , pEL )p*"»"Dl [2] 

RiS; ;EL)*PV p 


SIDE AFFECTED 
ELEMENTS AFFECTED 
PRSCRBD. VALUES AT INTC. PTS , EL 
MODIFY THE R MATRIX 


BF 

roi B+RZ BF S * C * RS 

Cl] p calculates’ THE FIELD VARIABLE CRADIANT INTERPOLATION MATRICIES 
C2] flS*l+pSC2; ;] 

[3] B+ 2 1 3 4$(tfE,pC)pOS[2; ;3 ,C1.53S[3; ;3 

[4] B+UNV R Z JACOB S)INPROD B 


CHKCNV 

CO] 

[ 1 ] 

[23 

[3] 

[4] 

[5] 

[63 

[73 

[83 

[93 

[103 

[113 

[123 

[133 

[143 


CNV+CHKCNV NIT 

p CHECKS CONVERGENCE OF EF , UP, VF 
p CNV RETURNS 1 IF ALL CONVERGED 


AND PF WITH LAST ITERATION VALUES 
0 IF NOT ALL CONVERGED 


CNV+ERR CHKCNV LSI 'EF' p 
+{FC=Q)/£ND fl 

CNV+CNVaERR CHKCNV LSI 'UF' p 
CNV+CNVaERR CHKCNV LSI 'VF' p 
CNV+CNV*ERR CHKCNV LSI 'PF' p 
END:*(CNV-Q)/itIMIT p 

“ STATUS ' CONVERGED WITH' NIT ' ITERATIONS ' p 

-►0 p 

LIMIT '.*(NIT<CIL)/0 p 

STATUS 'NOT CONVERGED IN' NIT ' ITERATIONS ' p 


CHECK EF 

JUMP IF NO FLOW CALCS. 
CHECK UF 
CHECK VF 
CHECK PF 

CHECK FOR CONVERGENCE 
MESSAGE TO USER 
EXIT 

EXIT FOR MORE ITERATIONS 
MESSAGE TO USER 


p 
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CHKCNVLS 1 

CO] CNV+ERR CHKCNV&S1 FLDNM ;AFC; FDD ; FLDSTAR ; IS 
Cl] o CHECKS FIELD VALUES FOR CONVERGENCE 

[2] p 


[3] CNV + 0 p 

[4] FLD**FLDNM p 

C5] FLDSTAR+lFLDNM , 'STAR' P 

C6] LS+(\FLD)<ERR*\ /\FLD p 

C7] LS+LSv(\PLD)<lE~10 p 

[8] LS+LSvl \ FLDSTAR )<ERRxf /\ FLDSTAR 

C9] LS*-LSv 0 = FLDSTAR p 

[10] +(*/LS)/QANT a 

[11] FLD*-(LS*-~LS)/FLD p 

[12] FLDSTAR+LS/ FLDSTAR p 

[13] AFC+\ (FLD-FLDSTAR)*FLDSTAR a 

[14] CNV+x/ERR>AFC a 

[15] *(CNV )/ 0 p 

[16] ♦ (ITER<CIL )/0 a 

[17] STATUS 'MAX ERROR '.FLDNM,': ', 

[18] +0 p 

[19] QANT; a 

[20] STATUS 'UNKNOWN CONVERGENCE OF' 

[21] CNV+ITERxl a 


DEFAULT IS NOT CONVERGED 
FIELD VALUES 

FIELD VALUES OF LAST ITERATION 
LOCATIONS WITH SMALL VALUES 
OR SMALL ABSOLUTE VALUES 
p OR RELATIVE SMALL STAR VALUES 

OR ZERO FLDSTAR VALUES 
CHECK IF NO VALUES WILL LEFT 
IGNORE SMALL FIELD VALUES 
AND THOSE SMALL STAR VALUES 
ABSOLUTE FRACTIONAL CHANCE 
RETURN A 1 IF ALL CONVERGED 
EXIT IF ALL CONVERGED 
EXIT IF NOT AT ITERATION LIMIT 
•T /AFCa DISPLAY MAXIMUM ERROR 
EXIT 

CAN'T DETERMINE CONVERGENCE 
FLDNMo WARN USER 

OK IF NOT FIRST ITERATION 


CHKEF 

[ 0 ] REC*-CHKEF ; RCD ; RGE 

[1] p CHECKS ENERCY FIELD VALUES TO SEE IF WITHIN PROPERTY DATA RANGE 

[2] p 

[3] RGD+MINMAX .PRQPDATA12 ; ;] p RANGE OF DATA VALUES 

[4] RCE^MINMAX EP a RANGE OP ENERGY FIELD VALUES 

[5] REC*-* /RGE< . iRCD p CHECK IF WITHIN RANGE 


CHKINPRES 

[0] CHKINPRES; A 

[1] p CHECK PRESRCRIBED PRESSURE NODE “PRNODE” SPECIFIED IN INPUT 

[2] p 

[3] +(2*CWC 'PRNODE ' )/0 p CHECK IF CONFLICT WITH MATL. PROP. 

[4] A+MINMAX, (PROPDATALSEL \ 6 ) [ 1 ; ;] p MIN. AND MAX. PRESSURE 

[5] *(.Al21<PRNODEi21 )/£RRl a CHECK IF TOO HIGH 

[6] +(il[l] >PRNODE [ 2 ] )/£RR2 p CHECK IF TOO LOW 

[7] -0 

[8] {Ml: 1 CHECK INPUT: TOO BICH PRESSURE SPECIFIED FOR PROPERTY DATA' 

[9] STOP 

CIO] £RR2:< CHECK INPUT : TOO LOW PRESSURE SPECIFIED FOR PROPERTY DATA' 

[11] STOP 
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CHKINPUT 

[0] CHKINPUT', A ;M 

[1] « CHECK INPUT PARAMETERS 

[l] " rHKTNPRES p CHECK PRESSURE SPECIFICATIONS 

4] f (~( 5<A ) vO . 2 >AW ( r /RZR )-L//?Z/?)/0 « CHECK ELEMENT ASPECT RATIOS 

[5] 'VAfWJiVC: ASPKCr tfAr/O 0F elements GREATER THAN 5' 


CfltfSS 

[ 0 ] 

[ 1 ] 

[ 2 ] 

[3] 

[4] 

[5] 
[63 
[7] 


OK+CHKSS CIC 
q CHECK IF STEADY STATE 


HAS BEEN REACHED 


OK* 0 p 

-*-(SSC=0)/0 p 

ok+*i\-cic p 

*(~OK)/U p 

STATUS ’STEADY STATE SOLUTION REACHED’ 


P 


INITIALIZE RETURN VARIABLE 
STEADY STATE CONTROL PAR AM . 
CHECK LAST ITERATION COUNTS 
EXIT IF ITERATIONS > 1 
MESSAGE TO USER 


CHKTF 

[0] RTC*CHKTF RGB ; RGT ; TMP ; TOL 

[1] p CHECKS TEMPERATURE FIELD VALUES 
C2] P 

[3] TOL*- . 1 p 

[4] RTC*1 p 

[5] TMP *- 1 FSTCM BC a 

C6] TMP*-\ FSTCM TMP a 

[7] TMP*-RN FSTCM TMP p 

[8] TMP *- 0 lvTMP p 

[9] -(2>l+prMP)/0 p 

CIO] RCB*-MINMAX*“9"TMP C : 1 ] o 
Cil] RGB*- (MEAN RGB)* 1 1* ( 1 +TOL )*0 . 5* 
C 12] RGT *-M INMAX TP a 

C13] RTC*-*/RCT< ,<,RGB p 


IF WITHIN PRESCRIBED VALUES 

TOLERANCE RANGE 
ASSUME OK 

ENERGY BOUNDARY CONDITIONS 
ONLY PRESCRIBED BCS. 

ONLY THIS REGION 
ANY SIDE 

EXIT IF NOT ENOUGH INFO. 
PRSCRB . TEMP . BCS . 

■/<t>RCB p ADD TOLERANCE TO RANGE 

RANGE OF TEMPERATURE VALUES 
CHECK IF WITHIN RANCE 


CLEANUP 

CO] CLEANUP',NAMES;ALTALP',C 

Cl] p ERASES VARIABLES NAMED IN QLV EXCEPT 

C2] p 

C 3 ] ->(0=0 NC ’QLV’)/ 0 p 

C4] NAMES * ( (1 +PQLV) ,10) +QLV p 

C5] ALTALP*’ABQDEFQHIiKLHNQPQSSIUYMMZ' " 

C6] NAMES* (~v7 v /NAMES 0 . -ALT ALP) /-NAMES p 

C7] C*DEX NAMES a 

C8] +Cl=x/C)/0 p 

C9] ’ VARIABLES NOT ERASED IN CLEANUP' p 

CIO] NAMESL ((~C)/(tpC )*C=0 ) ; ] p 


THOSE WITH ALT. CHARACTER NAMES 

EXIT IF QLV NOT AVAILABLE 

GET VARIABLE NAMES FROM QLV 

ALTERNATE ALPHABET 

DO NOT ERASE ALT. CHAR. NAMES 

ERASE THE REST 

CHECK IF ALL WERE ERASED 

MESSAGE TO USER 

DISPLAY THOSE NOT ERASED 
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CLS 

[ 0 ] 

[ 1 ] 

[ 2 ] 

[3] 

C4] 

[5] 

[ 6 ] 
[7] 
C83 
C9] 
[ 10 ] 
Cll] 
[ 12 ] 


CLS :RC;CTLS; DATS 
f> CLEARS 3270 SCREEN 
RC + 120 OSVO 2 4p 'CTLSDATS' 
+(v/2*RC)/N0SHARE 
RC+1 010 DSVC • CTLS ' 

CTLS+' PACE +1* 

RC*CTLS 

♦<0a.=«C)/O 

'RETURN CODE OP ' , (tRC). < FROM AP 120' 
DATS 

RC+QSVR 2 Up 'CTLSDATS' 

*0 

NOSHARE: 'OFFER TO AP 120 NOT ACCEPTED’ 


COLM 

[0] R+I COLM X 

[1] s RETURNS VALUES FOR INDEX 

[2] q 

[3] /?«-(J=i - ltpX)/X q 

[4] -*-(l<pp,I)/0 a 

[5] R+( UpRlpR a 


COORXE 

[0] COORXE: A 

[1] a CENERATES XI -ETA COORDINATES 

[2] q 

[3] il+ - l + Cl + il+2xJVC)x2 + 2xAr£> q 

[4] 4<-(2pp>! )pA a 

[5] , [0 . 5]<S$4 q 

[6] XEN*{2 ,0 .5** / pA)pA q 


OF LAST DIMENSION OF "X" 

SELECT DATA 

IF MULTIPLE COLUMNS, FINISHED 
ELSE RESHAPE, ELIM. LAST DIMENSION 


THE NODES 

NORMALIZED NODES LOCATIONS 
XI COORDINATES 
ADD ETA COORDINATES 
COORDS. WRT REGION COORDS. 


CPUCHK 


CO] 

CK+CPUCHK 


Cl] 

C2] 

q RETURNS A 0 IF CPU TIME LIMIT 
q 

IS EXCEEDED 

C3] 
C4] 
C 5] 

CK+~CPUTIMEZCPULIM q 
*(CK= l)/0 q 

STATUS 'CPU LIMIT EXCEEDED' q 

CHECK CPU TIME 
EXIT IP OK 
MESSAGE TO USER 

CPUTIME 


CO] 

CPU+CPUTIME 



[1] q RETURNS CPU SECONDS USED SINCE " TSTART •• WAS ISSUED 

[ 2 ] CP0^0.001xCWJ-[2]-ri[2] 
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DFN 

CO] 

[13 

[ 2 ] 

[3] 

[4] 
C5] 
[ 6 ] 


DF+DFN F;A;B;S 

o CALCULATES DERIVATIVES OF FIELD VARIABLE' AT THE ELEMENT NODES 


A +~ 1 "1 0 "1 1 “1 1 0 1 1 0 1 "1 1 
S+QLGSHPS 9 2p A A 
B+RZE BF S A 

Df+a 2 1*+/B*0MX<2 H)(.FIELN1) a 


1 0 0 0 a NODE XI -ETA COORDINATES 
SHAPE FUNCTIONS AT NODES 
GRAD. INTERP. MATRIX 
DERIVATIVES AT NODES 


DIAG 

CO] A-^DIAG X 

[1] a FORMS DIAGONAL MATRIX FROM A VECTOR X 

[2] 4*0 _ l + (-ipX)<K(2ppX)pO),X 


DIST 

[0] R*-A DIST B 

[13 a CALCULATES THE DISTANCE BETWEEN CARTESIAN POINTS A AND POINTS B 
[2] a LAST DIMENSION OF A AND B IS 2 COLUMNS OF X AND I'S 
[3 ] a 

[4] R+- (+/(S-A)*2)*0.5 a SQUARE ROOT OF SUM OF DIFF . SQUARED 


DX 

CO] 

[13 

[ 2 ] 

[33 

[4] 


R+DX T 

o RETURNS DERIVATIVES OF THE SHAPE FUNCTIONS WRT . X 
q T IS THE FIELD TYPE < E.C . 1- ENERGY) 

Q r+dDXTZFELTIT] 3 a SELECT FROM GLOBAL NESTED ARRAY 


DY 

[0] R+DY T 

[13 A RETURNS DERIVATIVES OF THE SHAPE FUNCTIONS WRT . Y 
[2] a T IS THE FIELD TYPE ( E.G . 1 -ENERGY) 

[33 A 

[4] R+oDYTIFELTZTn a SELECT FROM GLOBAL NESTED ARRAY 


EGYASN 

[03 LPV EGYASN F 

[13 a ASSIGN ENERGY FIELD VALUES 

[23 « 

[33 F+LPV\F A 

[4] FKPFEQZlil n+PFEQt2;l a 
[53 EF+FAN F A 


EXPAND FOR PRESCRIBED VALUES 
REINSERT PRESCRIBED VALUES 
ENERGY VALUES AT ALL NODES 
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EGYBAR 

[OJ EGYBARiA;B;CDT 

[1] a FORMS ENERGY EQUATIONS FOR SOLVING TRANSIENT RESPONSE 

[2] A 

[3] CEQ+LUMP CEQ a USE LUMPED CAPACITANCE 

[u] CDT+CEQ*TCTLL 3] a CAPACITANCE «■ LTIME STEP 

[5] KBAR+{THETA*KEQ)+CDT a NEW STIFFNESS MATRIX 

[6] A+(CDT+KEQ*THETA-1)+.*EPHI a INTERMEDIATE CALCULATION 

[73 B+, (RELAST*1-THETA)+THETA*EREQ a INTERMEDIATE CALCULATION 

[83 RBAR+A+B a NEW RESULTANT VECTOR 


EGYEVL 




[03 

EGYEVL ; KAE \ 

\CM\RC\RP 



[13 

a EVALUATES 

ENERGY MATRICES 



[23 

A 




[33 

PFEQ + 3 OpO 

0 

INITIALIZE 

PRESCRIBE VALUES MATRIX 

[43 

EGYMAT a 


CONSTRUCT , 

FLUID SUB-MATRICES 

[53 

FMFEQL 1 a 


FOAM FIELD 

EQUATIONS 


E GY MAT 

[03 EGYMAT 

[13 a CONSTRUCT THE ENERGY ELEMENT MATRICES 
[23 a 


[3] 

FLDMATLCM a 

FORM 

CM MATRIX 

[43 

FLDMATLRC a 

FORM 

RC VECTOR 

[53 

■* ( FC = 0 ) / 0 A 

EXIT 

IF NO FLOW CALCS 

[63 

FLDMAT LKAE a 

FORM 

KAE MATRIX 

[73 

MiNMrAflP a 

FORM 

RP VECTOR 


ENERGY 




[03 

ENERGY ; LPV ; FT ; RESULT 


[13 

a FORMS AND SOLVES 

THE 

TRANSIENT ENERGY EQUATION 

[23 

A 



[33 

+<EC=0)/0 a 


EXIT IF NO THERMAL ENERGY CALCS . 

[43 

FT+ 1 a 


FIELD TYPE IS ENERGY 

[53 

EGYEVL a 


FORM THE ENERGY EQUATIONS 

[63 

EGYBAR a 


MODIFY EQUATIONS FOR TRANSIENT FORM . 

[73 

PRSIDE FT A 


PRSCRIBE REGION SIDE BC'S 

[83 

PflSCflS a 


MODIFY FLUID ELEMENT EQS . 

[93 

LPV+PRSVLC FT a 


LOCATIONS OF PRESCRIBED VALUES 

[103 

LPV 

A 

SOLVE THE FLUID EQUATIONS 

[113 

LPV EGYASN RESULT 

A 

ASSICN THE FIELD VALUES 



135 


FAN 

[ 0 ] 

Cl] 

[ 2 ] 

[3] 

[4] 
C 5 1 
[ 6 ] 

[7] 

[ 8 ] 


R+PAN F’.MAX 

fl RETURNS field variable at all nodes 

a F IS FIELD VARIABLE 

n R IS MAXIMUM OF NNG FOR ALL FIELD TYPES 


R+-F a 

MAX+r /NNCxpFELT p 
■* (MAX-pF) /0 p 
R+LINFV F p 


INITIALIZE RETURN VARIABLE 

MAXIMUM NODES PER ELEMENT 

EXIT IF ALREADY SAME AS LARGEST ELEMENT 

CONVERT LINEAR TO QUADRATIC VALUES 


FLDMATLCM 

[0] FLDMATLCM-, A 

[1] p FORMS CM MATRIX FOR PL DM AT 

[2] p 

[3] A«-l FVEB RHO 

[4] A«-l 2 4 3^ (NS 1)*0PAX(2 3 )A 

C 53 A*- (NS DINPROD A 

[6] CM+NODEM AINTGRT A 


PLDMATLKA 

[0] FLDMAT&KA-.A 

Cl] p FORMS KA MATRIX FOR FLDMAT 

[2] P 

[3] A*(NS 2 )*0PAX ( 2 3 ) C ( 2 FVEB RH0)*2 FVEB UF) 

[ 4 ] KA+NODEM AINTGRT A IN PROD 1 2 4 3 6)£)X 2 

[5] A+(NS 2)*OPAX(2 3)((2 FVEB RHO)* 2 FVEB VF ) 

[6] KA+KA+NODEM AINTGRT A INPROD 1 2 4 3S)0Y 2 


FLDMATLKAE 

[0] FLDMAT LKAE-, A 

[1] p FORMS KAE MATRIX FOR FLDMAT 

[2] p 

[3] A*-(NS 1)*0PAX(2 3 ) ( ( 1 FVEB RH0)*1 FVEB UF) 

[4] KAE+NODEM AINTGRT A INPROD 1 2 4 3 $£>X 1 

[5] A*(NS l)*OPAX(2 3 ) ( ( 1 FVEB RHO )* 1 FVEB VF) 

[6] KAE+KAE+NODEM AINTGRT A INPROD 124 3$DY 1 
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FLDMATEKNN 
CO] FLDMATEKNN ; 4 

[1] o FORMS Kll K 12 X21 K22 MATRICIES FOR FLDMAT 
C2] o 

[3] 4*1 2 4 3$ (DX 2)*0P4X(2 3)(2 FVEB VIS ) 

C4] Kll+NODEM AINTCRT (DX 2)INPR0D A 
CS] K12+N0DEM AINTCRT (DY 2 )INPROD A 

[6] 4*1 2 4 3 $(£>X 2)*0P4X(2 3) (2 FVEB VIS ) 

[7] K21+N0DEM AINTCRT (DX 2)INPR0D A 

[8] K22+NODEM AINTCRTiDY 2)INPR0D A 


FLDMATELCN 

[0] FLDMATELCN-, A 

Cl] o FORM LC\ AND LC2 MATRICIES FOR FCDMAT 
C2] p 

C3] A*(.NS 4 )INPR0D 1 2 4 3 WX 2 
C 4 ] LC1+NODEM AINTCRT 4*0P4X(2 4) (2 FVEB RHO ) 
C 5] 4*(tfS 4 )INPR0D 1 2 4 3$py 3 
C6] LC2*-N0DEM AINTCRT 4*0P4X(2 4)(3 FVEB RHO ) 


FLDMAT ELN 
CO] FLDMAT ELN 

Cl] p FORMS LI AND L 2 MATRICIES FOR FLDMAT 
C2] p 

C3] L\*-N0DEM AINTCRTiDX 2)INPR0D 1 2 4 3 IstNS 4 

C4] L2+-NODEM AINTCRKDY 3)INPR0D 1 2 4 3$JVS 4 


FLDMAT EMM 
CO] FLDMAT EMM ; 4 

Cl] p FORMS MASS MATRIX FOR FLDMAT 
C2] p 

C3] 4*2 FVEB RHO p 

C4] 4*(/VS 2 )INPR0D 1 2 4 3$(JVS 2)*0P4X(2 3 )4 p 

C5] MM+NODEM AINTCRT A a 


FLDMAT ER 

CO] FLDMAT ER\GV 

Cl] p FORM R VECTORS FOR FLDMAT 

C2] p 

C3] CV+GRAVVEC 

C4] RU+CVll ; ; 1+NODEV SINTCRT BDY 2 
C5] RV+CV 12 \ ;1*N0DEV SINTCRT BDY 3 


RHO AT QUAD. ELMNT NODES 
INTER. CALC. 

MASS ( INERTIA ) MATRIX 
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FLDMATLRC 

[0] PLDMAT LRC 

[1] a FORMS RC VECTOR FOR PLDMAT 

[23 « 

[3] RC*-NODEV SINTCRT BDY 1 

[4] RC+RC+RCOND 


FLDMATLRP 
[03 FLDMATLRP 

[13 p FORMS AND ASSIGNS RP MATRICIES 
[23 « 

[33 RP*-RPRES 


FOR PLDMAT 


PLDSLV 

[03 F+FLDSLV LPV-.X 

[13 n SOLVES EQUATIONS AND ASSIGNS 

[23 a 

[33 *(.ICTL = 0 1 )/SS,NR « 

[43 £S: p 

[53 F+.RBARftKBAR a 

[63 -*-0 p 

[73 NR: p 

[83 X+LPV / FLDSLV C.STAR A 

[93 F+(KBAR+ .*X)-RBAR p 

[103 F*-X+(.-F)MBAR p 


FIELD VALUES 

CHECK TYPE OF ITERATION METHOD 
SUCCESSIVE SUBSTITUTION 
SOLVE SIMULTANEOUS EQUATIONS 

EXIT 

NEWTON -RAPHSON 

FIELDS FROM LAST ITERATION 

NEW FORCING FUNCTION 

SOLVE SIMULTANEOUS EQUATIONS 


FLDSLVLSTAR 
[03 F+FLDSLVLSTAR 

[13 p RETURNS FIELDS FROM LAST ITERATION 
[23 p 

[33 *(~FC=0)/W1 p 

[4] F+EF o 

[5] +0 a 

C6D Nl:EF,UF,VF,PF n 


CHECK IF PLOW CALCS. WERE SOLVED 
FOR CASE OF NO FLOW CALCULATIONS 
EXIT 

FOR CASE WITH PLOW CALCULATIONS 


138 


FLOW 

CO] FLOW iLPV; FT; RESULT 

Cl] o FORM AND SOLVE THE FLUID FLOW 
C2] A 

C3] +(FC- 0)/0 A 

[4] FT* 2 3 4 a 

C5] FLWEVL a 

C6] FLWBAR A 

C7] PRSIDE FT a 

[8] PRSPRES A 

C9] PRSCRB a 

CIO] LPV*PRSVLC FT a 

Cll] RESULT+PLDSLV LPV a 

C12] LPV FLWASN RESULT a 


EQUATIONS 

EXIT IF NO FLOW CALCS. 

FIELD TYPES U. V, AND P 

FORM THE FLUID EQUATIONS 

MODIFY EQUATIONS FOR TRANSIENT FORM. 

PRSCRIBE REGION SIDE BC'S 

PRESCRIBE PRESSURE NODE 

MODIFY FLUID ELEMENT EQS. 

LOCATIONS OF PRESCRIBED VALUES 
SOLVE THE FLUID EQUATIONS 
ASSIGN THE FIELD VALUES 


FLWASN 

CO] LPV FLWASN P 

Cl] a ASSIGN FLUID FLOW FIELD VALUES 
C2] A 


C3] F*LPV\F a 

C 4 ] Fl(.PFEQ:i;l )l+PFEQZ2il a 

C5] UF+ (A+NNG 2)+F a 

C6] F+A+P a 

C7] VF+{A+NNC 3 ) +F a 

C8] F*A*F a 

C9] PF* ( NNG 4 )+F a 

CIO] (UF VP PP)*FAN"UF VP PF a 


EXPAND FOR PRESCRIBED VALUES 
REINSERT PRESCRIBED VALUES 
EXTRACT U VELOCITY VALUES 
DROP THOSE VALUES 
EXTRACT V VELOCITY VALUES 
DROP THOSE VALUES 
EXTRACT PRESSURE VALUES 
EXPAND TO ALL NODES 


FLWBAR 

CO] FLWBAR; A ;B;CDT 

Cl] a FORMS FLOW EQUATIONS FOR SOLVING 
C2] a 

C3] CEQ+LUMP CEQ a 

[4] CDT+CEQ*TCTL13] a 

C 5] KBAR* ( THETAxKEQ ) +CDT a 

[6] A*(CDT+KEQ*THETA-1 )+. xPPHI a 

C 7 ] B* , ( RFLASTx 1 - THETA ) +THETA *FREQ a 

C8] RBAR+A+B a 


TRANSIENT RESPONSE 

USE LUMPED CAPACITANCE 
CAPACITANCE ♦ AT IMS STEP 
NEW STIFFNESS MATRIX 
INTERMEDIATE CALCULATION 
INTERMEDIATE CALCULATION 
NEW RESULTANT VECTOR 


FLWEVL 

CO] FLk/EV LiMMiKA;KC;Kll;K12;K21;K22;Ll;L2iRC:RU:RV:LCl'LC2'RP 
* EVALUATES FLUID FLOW ELEMENT MATRICES '“ V ' LC1 ' LC2 ' RP 

P 

FLWMAT ° P0 " INITIALIZE PRESCRIBE VALUES MATRIX 

pmpfoai \ CONSTRUCT FLUID SUB-MATRICES 

FMFEQal P FQRH FIELD EQUATIONS 


Cl] 

C2] 

C3] 

C4] 
C 5] 
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FLWMAT 

CO] 

FLWMAT 

[1] 

a CONSTRUCT Tl 

[2] 

a 

[3] 

FLDMATLMM a 

C4] 

FLDMATLKA a 

[5] 

FLDMATLKNN a 

[6] 

FLDMATLLN a 

[7] 

FLDMATLLCN a 

[8] 

FLDMATLR a 

FMFEQL 1 

[0] 

FMFEQbl 

Cl] 

a FORM ENERGY 

[2] 

Q 

C 3 ] 

CEQ+CM a 

[4] 

EREQ+RC a 

[5] 

-(5(7=1 )/FL0W 

[6] 

KEQ+ (pCEQ )p 0 

[7] 

-►0 a 

[8] 

FLOW:KEQ*KAE i 

[9] 

EREQ+EREQ+RP 


FORM MASS MATRIX 
FORM KA MATRIX 

FORM Kl\ K12 K2i K22 MATRICIES 
FORM LI L2 MATRICIES 
FORM LC1 LC2 MATRICIES 
FORM RU RV VECTORS 


CAPACITANCE MATRIX 
ENERGY RESULTANT VECTOR 
CHECK FOR FLUID FLOW CALCS. 
NO FLOW STIFFNESS MATRIX 
EXIT 

WITH FLOW STIFFNESS MATRIX 
WITH FLOW RESULTANT MATRIX 


FMFEQL 2 

[0] FMFEQL2 ; Z1 ; Z2 ; Z3 

[13 a FORM FIELD EQUATIONS 

[2] a 

[3] Z1 *(NNG 2 4)p0 o 

[4] Z2*- (NNG 2 2 )p0 a 

[5] Z3*-(NNG 4 4 ) p 0 a 

[63 KEG<-(K4+K22 + (4 + 3)xK 11), X12.il a 

[73 KEQ+KEQ, [13X21 , (XA+Xll+(4*3 )*K22 ) ,L2 

[83 KEQ+KEQ, [13LC1 ,LC2 ,Z3 a 

[93 FREQ*RU .ill RV .ill ((NNG 4),l)p0 a 

[103 CEQ*-M M.Z2.Z1 a 

[113 CEQ+CEQ.i 11Z2.MM.Z1 a 

[123 CEQ+CEQ , [13 ($Z1), (#>Z1),Z3 a 


ZERO MATRIX 
ZERO MATRIX 
ZERO MATRIX 

X MOMENTUM ( STIFFNESS ) 
I MOMENTUM ( STIFFNESS ) 
CONTINUITY ( STIFFNESS ) 
FLOW RESULTANT VECTOR 
X MOMENTUM (INERTIAL) 

Y MOMENTUM ( INERTIAL ) 
CONTINUITY ( INERTIAL ) 
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FRONTEND 

[0] FRONTEND A 

[1] p PERFORMS UPFRONT 

[2] P 

[3] TSTART P 

[4] CLEANUP p 

[5] STATMSG * 0 Op' ' p 

C6] ND+A p 

[7] INPUT P 

C8] CRIDGEN p 

C9] COORXE p 

CIO] MAP a 

Cll] INITIAL a 


ONCE ONLY FUNCTIONS 

START CLOCK FOR TIME CHECKING 

ERASE OLD VARIABLES 

INITIALIZE STATUS MESSAGE 

NUMBER OF DIVISIONS PER SIDE OF REGION 

SPECIFY INPUT PARAMETERS 

GENERATE ELEMENT NODE NUMBERING 

GENERATE XI-ETA COORDINATES OF THE NODES 

MAP XI-ETA NODE COORDINATES INTO X-Y SYSTEM 

INITIALIZE VARIABLES 


FSTCM 

CO] R+X FSTCM M 

Cl] p RETURNS SUB-MATRIX WHERE MEMBERS OF X MATCH FIRST COLUMN OF M 
[2] p 

C3] R* (M[ ; 1 ] eX )/0 1+M 


FST 1 

CO] R*FST 1 A 

Cl] p PUTS 1 IN COLUMN OP EACH ROW OF A WHERE FIRST POSITIVE CHANGE OCCURS 
L 2 J R+< \A 


FVEB 

CO] 

[ 1 ] 

[23 

[3] 

C4] 

C53 

[ 6 ] 

[7] 

C8] 

C9] 


R<-FT FVEB X ; ET 

p RETURNS THE FIELD VARIABLE ON AN ELEMENT NODE BASIS 
p X IS FIELD VAR . ON A GLOBAL NODE BASIS FOR HIGHEST ORDER ELEMENT 
fl FT IS THE TYPE OF FIELD ( E.C . 1 -ENERGY ) 

p 


ET+FELTLFTl P 
*(,ET-1 2'i/TYPE\,TYPE2 p 
TYPEUR+XtELNliX 3 5 7]] 

-*■0 p 

TYPE2:R+XtELNl p 


CONVERT FIELD TYPE TO ELEMENT TYPE 
CHECK TYPE OF ELEMENT 

FIELD VARIABLES ON LINEAR ELEMENT BASIS 
EXIT 

FIELD VARIABLES ON QUAD. ELEMENT BASIS 
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FVGB 

CO] 

[1] 

[ 2 ] 

C3] 

[4] 

[5] 

[6] 
[7] 
C8] 

[ 9 ] 

[ 10 ] 
[ 11 ] 
[ 12 ] 


R+FT FVGB X\A\ET 

a RETURNS THE FIELD VARIABLE ON AN GLOBAL NODE BASIS 
n X IS FIELD VAR. ON A GLOBAL NODE BASIS FOR BICHEST ORDER ELEMENT 
n FT IS THE TYPE OF FIELD ( E.G . 1 -ENERGY) 


ET+FELTZFT1 a 
*(ET-1 2)/2'1.2’2 n 
£1: a 

X+(2pA+l+2*ND)pX a 
X«-(A<-Apl 0)/X a 
R*- ,A/X a 
-»0 a 

T2:R+X a 


CHANCE FIELD TYPE TO ELEMENT TYPE 

CHECK TYPE OF ELEMENT 

CONVERT TO LINEAR GLOBAL BASIS 

MAKE FIELD VARIABLE A MATRIX 

REMOVE LINEAR NODE COLUMNS 

REMOVE LINEAR NODE ROWS 

EXIT 

LEAVE ON QUADRATIC GLOBAL BASIS 


FVIP 

[0] R+FVIP F 

[1] a CALC FIELD VARIABLES AT INTEGRATION POINTS 
C2] a 

[3] /?■*■ + / (NS 2 ) [ ; $ ; 13 *0PAX(.2 3)FZELN1 


CRAW EC 

[0] CV+CRAVVEC',R;RZ 

[13 a RETURNS GRAVITATIONAL BODY FORCE VECTOR 


[2] A 

[3] R+RHO a 

[4] RZ*RH0ZER0 a 

[5] RZIEA1 SFe 1 2 1+RZ a 

[6] R+R*1-R*RZ a 

[7] GV*- ( 2 FVEB R)*0PAX( 2 3 ) (.NS 2) o 

[8] CV+-NODEV AINTCRT CV a 

[93 GV+(.GRZZn*CV),Z0.51CRZZ21xCV a 


DENSITY 

REFERENCE DENSITY 

SET SOLID AND 2$ NODES TO RHOZERO 
NORMALIZED DENSITY 
MULTIPLY BY THE SHAPE FUNCTIONS 
INTEGRATE ON GLOBAL BASIS 
MULTIPLY BY GRAVITY 


CRIDAELN 


[ 0 ] 

C13 

[23 

[33 

[4] 

[53 

[63 

[73 

[83 

[93 

[103 


ELN+CRIDbELN CNQ-.A 

a GENERATES QUADRATIC NODE NUMBERS ( ELEMENT BASIS ) FOR CRIDCEN 

A 


A<-1 , ( - l + 3*A7D)pl 10 a 
ELN*A\CNQ A 
ELN+-A \ELN A 

ELNZ:Al+ELNZl~l*A+IEAl~Al a 
ELtf[A;3+FLff[ - l+A;3 a 
ELN + 1 3 2 4$(4piW>,3 )pELN A 
ELN+.Z 1 23 , [3 43 ELN a 
ELN*ELNZ ; 7 8963214 53 


EXPANSION VECTOR 

EXPAND CLOBAL NODE NUMBER ON COLUMNS 

EXPAND GLOBAL NODE NUMBERS ON ROWS 

DUPLICATE COLUMNS 

DUPLICATE ROWS 

RESHAPE NODE NUMBERS 

AND SEPARATE FOR EACH ELEMENT 

REORDER FOR STANDARD ELEMENT 
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CRIDALNODE 

[0] LNODE+GRIDALNODE GNL ;A;N 

[1] fl GENERATES LINEAR NODE NUMBERS ( ELEMENT BASIS ) FOR CRIDGEN 
C2] fl 

[3] A+2*~19,*(2,N')piN+NDri 

[4] A+<.((N,N)piN*2)l;Al )CA;] 

[5] A4-((NDx2),ND,2)p((.2xNE) ,2)pA 

[6] A+( (ND*2 ) ,4 )pl 3 2*(W>,4,tfO)pl 3 2M 
C7] LNODE*A C : 3 421] 


GRID ARSE 

C 0 ] RSE+GRIDARSE ;A-.EN 

Cl] n GENERATES REGION SIDE ELEMENT 

[2] B 

[3] EN+(2pND)p\ND*2 b 

[4] i5+EAf[l f ltp£'JV;] b 

[5] 4+Jl,Cl]W/irCjl,"l*pffW fl 

[6] RSE*A C2 4 1 3;] b 


NUMBERS 

ELEMENT NUMBERS 
ELEMENTS ON SIDES 3 AND 1 
ADD ELEMENTS ON SIDES 4 AND 2 
ELEMENT NUMBERS- -EACH SIDE OF REGION 


CRIDARSLN 

[0] RSLN+CRIDARSLN CNL;A 

Cl] fl GENERATES REGION SIDE LINEAR NODE NUMBERS FOR CRIDGEN 
C2] b ' 

C3] i4«-2 + [l]"leCJV£ b LINEAR NODE NUMBERS FOR SIDES 1 AND 3 

C4] A*-A, Cl]$2+C2]"l<t>CiV£ b ADD NODE NUMBERS FOR SIDES 2 AND 4 

C5] RSLN+A C 1 3 2 4;] a LINEAR NODE NUMBERS ON EACH SIDE 


CRIDARSN 

CO] RSN+CRIDARSN CNQ-.A 

Cl] fl GENERATES RECION SIDE QUADRATIC NODE NUMBERS FOR CRIDGEN 
C 2] B 

C3] il-*-2t Cl] ~19GNQ a NODE NUMBERS FOR SIDES 1 AND 3 

C4] A+A , Cl]«t2 + C2]“l<t>CAf$ b ADD NODE NUMBERS FOR SIDES 2 AND 4 

C5] RSN+A tl 3 2 4;] b NODE NUMBERS ON EACH SIDE OF RECION 


NODE NUMBERING 


CRIDGEN 

CO] CRIDGEN', A \ GNL ; GNQ 

Cl] fl GENERATES ELEMENT 
C2] b 

C3] CNQ*(2pA)p\ (A*-l + 2xND)*2 fl 

C4] GNL*{2pA )p i (A+1+ND)*2 fl 

C5] RSN+CRIDARSN GNQ n 

C6] RSELMT+-GRI DARSE n 

C7] RSLN*GRIDARSLN GNL a 

C8] ELN*-CRIDAELN GNQ a 

C9] LNODE+CRIDALNODE GNL fl 


GLOBAL NODE NUMBERS- -QUAD . ELEMENTS 
GLOBAL NODE NUMBERS- -LINEAR ELEMENTS 
NODE NUMBERS ON EACH SIDE OF RECION 
ELEMENT NUMBERS --EACH SIDE OF REGION 
LINEAR NODE NBRS.-EACH SIDE OF RECION 
QUADRATIC NODE NUMBERS- -ELEMENT BASIS 
LINEAR NODE NUMBERS- -ELEMENT BASIS 
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HERON 

[0] ACHERON TiS 

[1] a CALCULATES AREAS OF TRIANGLES "T" USING HERON'S FORMULA 

[2] b LAST DIMENSION OF T CONTAINS 3 SIDE LENGTHS 

[3] S«-0.5x+/r 

[ 4 ] 4+(0rS*(S-l COLM n*(S-2 COLM T)* (S-3 COLM D)*0.5 
IEA1 

[0] R+IEA 1 A 

Cl] b RETURNS INDICIES OF LOCATION OF ALL ONES IN VECTOR A 
[2] R+R/ (,\pA)*R->-A = l 


IFST1 

COl R+IFSTX A 

[1] b RETURNS INDICIES OP FIRST POSITIVE CHANGE IN EACH RON OF A 
C2] R+(, A)/ ,A*(pA)p\~l + pA->-FSTl A.l 


INBDY 

[0] INBDY 

[1] fl SPECIFY RECION BOUNDARY CONDITIONS 

[2] b FIELD 

[3] b 1 TEMPERATURE 

[4] b 2 FLUID VELOCITY IN X DIRECTION 

[5] b 3 FLUID VELOCITY IN Y DIRECTION 

[6] b TYPES 

[7] o 1 PRESCRIBED VALUE 

C8] b 2 PRESCRIBED FLUX(SICN CONVENTION FLUX IN IS POSITIVE ) 

[9] BO 5 OpO 

[10] b FIELD TYPE REGION SIDE PROPERTIES 

[11] BOBC, 1 1 1 1 “.5 

[12] BC*-BC . 1 113 1.5 

[13] B 

[14] BOBC.2 114 0 

[15] BOBC, 2 1120 

[16] BOBC.2 1110 

[17] BOBC.2 113 0 

[18] B 

[19] BOBC.3 114 0 

[20] BOBC.3 112 0 

[21] BOBC.3 1110 

[22] BOBC, 3 1130 

[23] BO$BC 

[24] fi 

[25] PRN0DE+2 0 A PRESCRIBED PRESSURE NODE 


144 


INCEOH 

CO] INCEOH; A ;RR;ZR 

Cl] a SPECIFY GEOMETRY INPUT PARAMETERS 
C2] a 

[3] INREC a 

[4] CF + 1 a 

[5] CRZ*~2E~2 0 A 


INITIAL 

CO] INITIAL; A 

Cl] a INITIALIZE VARIABLES 

C2] A 

C 3 ] CRKINPUT a 

C 4 ] BN* 1 a 

C5] AINTCON AN 1*2 a 

C 6 ] SINTCON SNI* 2 a 

C7] AREA a 

C8] TIME*TCTL Cl] a 

C9] INITTEMP a 

CIO] INITVEL a 

C 1 1 ] INITPRES a 

Cl 2] INITMPROP a 

C 13 ] EREQ*RPRES+RCOND a 

C 14 ] 4-«-, Cl 2 1CRAVVEC a 

CIS] FREQ*A,ilU(NNC 4),l) p o a 

C16] SETLAST a 

C17] SETSTAR a 

C 1 8] SAVEFLDS*0p0 a 


SPECIFY COORDINATES OF REGION 
CARTESIAN COORDINATES 
GRAVITATIONAL CONSTANTS 


CHECK INPUT PARAMETERS 
REGION NUMBER 

AREA INTEGRATION CONSTANTS 
SURFACE INTEGRATION CONSTANTS 
VOLUME AND SIDE AREAS OF ELEMENTS 
INITIALIZE TIME 
INITIALIZE TEMPERATURE 
INITIALIZE VELOCITIES 
INITIALIZE PRESSURE 
INITIALIZE MATERIAL PROPERTIES 
ENERGY RESULTANT VECTOR 
GRAVITATIONAL BODY FORCE VECTORS 
FLOW RESULTANT VECTOR 
SET LAST TIME ITER. VARS. 

SET LAST CONVERGENCE ITER. VARS. 
INITIALIZE STORAGE OF FIELD VARS. 


INITMP 

CO] INITMP;A;DATA 

Cl] a INITIALIZE MATERIAL PROPERTIES : 
C2] a 

C3] TP* (NNG 2 )oTINIT a 

C4] DATA*(PROPDATALSEL 1 4 6) Cl 4 2; 

C5] A*(PF,ll.SlTF)INTQR DATA a 

C6] EF<-4C1;1:] a 

C7] UPPROP A 


CASE WHERE PRESSURE NODE PRESCRIBED 

INITIALIZE TEMPERATURES TO UNIT 
;] a ASSUME INITIALLY 1$ 

INTERPOLATE 

ENERGY FIELD 
UPDATE ALL PROPERTIES 


INITMPROP 
CO] INITMPROP 

Cl] a INITIALIZE MATERIAL PROPERTIES 
C2] A 

C3] +(2 =DNC 'PRNODE' )/NP a 

C4] STOP a 
C5] NP; INITMP a 


USINC AVERAGE CONDITIONS 

CHECK IF PRESSURE NODE PRESCRIBED 
STOP EXECUTION 
INITIALIZE PROPERTIES 
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INITPRES 

[0] INITPRES 

Cl] ft INITIALIZE PRESSURE FIELD 
C2] ft 

C3] *(0=Q NC 'PFINIT' )/NXTl ft 

C4] PF+PFINIT ft 
C5] >0 ft 

C6] NXT1:PP+(NNG DpPRNODE 11 ] ft 

C7] INITMP ft 

C8] PF+FAN STATPRES ft 


CHECK IF PFINIT EXISTS 
USE VALUES IN PFINIT 
EXIT 

FIRST USE REFERENCE PRESSURE 

CALC. DENSITY 

CALC. STATIC PRESSURES 


INITTEMP 

CO] INITTEMP i A 

[1] ft INITIALIZE TEMPERATURE 

C23 ft 

C 3 ] TINIT+,e$(2pl+A)pU+2xND)INTERVALS ’.5 1.5 ft INITIALIZE TEMPERATURE 


INITVEL 

[0] INITVEL 

[1] o INITIALIZE FLUID VELOCITIES 

[2] ft 

[3] -►(~a/2=G NC 2 6p * UFINITVFINIT' )/ZER0 ft CHECK FOR INITIAL VELOCITIES 

C4] UF+UFINIT o INITIAL U VELOCITY 

C 5] VF+VFINIT o INITIAL V VELOCITY 

C6] +0 ft EXIT 

[7] ZERO : UF+(NNG 2)p0 ft ASSUME ZERO U VELOCITY 

[8] VF+(NNG 2 ) p 0 ft ASSUME ZERO V VELOCITY 


146 


INPLYGN 
[ 0 ] 

[13 
C23 
[33 
[43 
[53 
[63 
[73 
[83 
[93 
[103 
[113 
[123 
[133 

[ 14 ] 

[ 15 ] 

[ 16 ] 

[ 17 ] 

[ 18 ] 

[ 19 ] 

[ 20 ] 
[ 21 ] 
[ 22 ] 

[ 23 ] 

[ 24 ] 


R+XY INPLYGN Q ;AREAP;CNPi INP-.MXY ;NXY ;NN;S;S1 ;S2 
FINDS IF POINTS XY ARE IN POLYGONS Q 

XY IS 2 COLUMN MATRIX AND Q p ( NBR . POLYGONS ) (AW?. NODES EA.) 2 
Q NODES ARE NUMBERED SEQUENTIALLY AROUND PERIMETER 
R RETURNS NBR. OF THE POLYGON CLOSEST TO EACH XY p (1+pXY) 
APPROACH IS TO FIND AREAS OF POLYGONS BY DIVIDING INTO TRIANGLES AND 
COMPARING TO SUM OF AREAS OF TRIANGLES OF XY AND EACH POLYGON SIDE 


AW«-(p«)[23 fl 

INP+INP , [ 1 . 5 ] 1 <t>INP+ \ NN a 

CNP*-Q [ ; INP ; ] n 

MXY+MEAN 1 3 2«)Q o 

MXY*2 1 3<*(.NN,pMXY)pMXY a 

S+CNPh ;1;)DIST CNPii ; 2 ; ] a 

Sl+CNPl ; ; 1 ; )DIST MXY fl 

S2*CNPL; ;2-,']DIST MXY a 

AREAP+ I +/HERON S,S1,[2.53S2 a 

CNP+Ul + pXY) ,pCNP)pCNP a 

JVXY-*-(l<»i4)S)(l<t>pCVP[; ; ; 1 ; ] )p$XY fl 

S«-((l+pXY),pS)pS b 

S1*-CNPC ; ; ; 1 ; JDIST NXY a 

S2+CNPI ; ; ; 2 ; ] DIST NXY a 

R+$\ +/HERON S,S1,[3.5]S2 fl 

R*(. (>HR)-(UipR)pAREAP)+OPAX 2 AREAP b 

R*IFST1U/R) = 0PAX 1 R a 


NUMBER OF NODES PER POLYGON 
INDS. OF SUCCESSIVE NODE PAIRS 
COORDS. OF SUCC. NODE PAIRS 
MEAN X Y COORDS. EACH PLYCN . 
RESHAPED FOR NBR. OF NODES 
SIDE LENGTHS OF POLYGONS 
LENGTHS OF 1 NODE TO MEAN PT . 
LENGTHS OF 2 NODE TO MEAN PT . 
AREA OF EACH POLYCON 
RESHAPE NODE PAIRS FOR EACH XY 
RESHAPE XY ALL PLYGNS . , NODES 
RESHAPE PLYCN SIDE LENGTHS 
CALC LENGTHS OF 1 NODE TO XY 
CALC LENGTHS OP 2 NODE TO XY 
AREAS EACH PLYCN. WITH EA XY 
NORMALIZED DIFF. IN AREAS 
CLOSEST PLYCN. FOR EACH XY 


INPROD 

[0] D+A INPROD B;C;RA;N-,CB;R;E 

[13 fl CALCULATES INNER PRODUCT OF MULTIPLE MATRICIES 

[2] fl 

[33 -*-((ppiO*(pp0))/CHir 

[4] -»(~(x/(ff-e2 + - 2epA )=(2 + “2ep0) ) )/CHK 

[5] -<•( (C«-l+"lep,A )*(R*-\*~29pB))/NXT 

[63 CHK: ' INCORRECT p IN INPROD' 

[ 7 ] STOP 

[8] NXT-.RA+\ + ~2$pA 

[9] Cfl«-l + “l$p0 

[103 D+(E,RA,CB)p+/L\) (4 2 3 Vi(.CB,N,RA ,C)pA )*3 2 1 **)(RA . (.N+*/E) .R.CB )pB 
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INPROG 

[0] 

INPROG 

[1] 

p SPECIFY PROGRAM C 

[2] 

p 

[3] 

CIL+ 6 P 

[4] 

CPULIM+ 90*3600 P 

C 5] 

TIMELIM+ 72x3600 P 

C6] 

EC+1 R 

C7] 

E/?J?«-0.05 p 

C 8 ] 

FC + 1 p 

[9] 

FELT+2 2 2 1 P 

[10] 

ICTL+0 p 

[11] 

SFN + 5 p 

[12] 

TCTL + 0 200 2 p 

[13] 

THETA+0.5 P 

[14] 

rsoo p 

[15] 

SSC+l p 

INPUT 

[0] 

INPUT 

[1] 

p SPECIFY AND PRINT 

[2] 

p 

[3] 

INGEOM R 

[4] 

INBDY P 

[5] 

INPROC R 

[6] 

NE+ND*ND R 

INREG 

[0] 

INREG-.RR-.ZR'.A 

[1] 

a SPECIFY INITIAL R 

[2] 

n 

[3] 

a' LARGE SQUARE BOX’ 

[4] 

RR+.1*0 0.5 1 1 1 

[5] 

ZR*- .1*0 0 0 0.5 1 

[6] 

RZR+-RR, [0. 5]Zfl 

INTERVALS 

[0] 

I*N INTERVALS P:B; 

[1] 

a GENERATES N INTER 

[2] 

fl-H+P 

[33 

E-*-l+P 

[4] 

I+B+(((B-B)*N)x~ 1 + 


CONVERGENCE ITERATION LIMIT 
LIMIT FOR CPU TIME 
TIME LIMIT 

ENERGY CALC. CONTROL (0 -NO THERMAL CALC. ) 
ALLOWABLE FIELD CONVERGENCE ERROR 
FLUID FLOW CALC. CONTROL (0 -NO FLOW CALC.) 
ELEMENT TYPES FOR EACH FIELD 
ITERATION CONTROL ( 0-SUB \-NEWTON-RAPHSON) 
SAVE FIELDS AT EVERY SFN TIME STEPS 
TIME START , END AND INCREMENT CONTROL 
TRANSIENT ALGORITHM CONTROL PARAMETER 
TIME STEP CONTROL (0-CONSTANT . 1 -VARIABLE) 
STEADY STATE CONTROL 


CEOMETRY FACTORS 
BOUNDARY CONDITIONS 
PROCRAM OPERATION PARAMETERS 
NUMBER OF ELEMENTS PER REGION 


Z COORDINATES OF REGION 
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INTQR 

CO] 

Cl] 

C2] 

C3] 

C4] 

C 5] 

C6] 

C7] 

C8] 


R+XY INTQR D;A;B;MAX;MIN;ND;NXY 



fl D ARE INTERP . DATA X , Y AND Z ( MAY BE MULTIPLE Z) p £3 NR 8 
« NR IS THE NUMBER OF QUADRATIC REGIONS 

b R RETURNS INTERP. VALUES WITH 2 PARTIAL DERIV . pC2 + l + p D) 3 (ItpXY) 


-►(0 = 1 +pXY)/N0XY B EXIT IP NO POINTS 

B 


[9] (4 MAX MIN)+N0RMR *. [2 3]D p 

[10] ND-*(pD )p(84 p 

[11] NXY-* (XY -OPAX 2(2+MIN))*OPAX 2(2*MAX-MIN) 

[12] 4 -*-3 1 2$4-<-/V0[l 2;;] a 

[13] A-*NXY INPLYCN A a 

[14] R+ND INTQRLREC NXY ,4 p 

[15] R+INTQRbUNN (R MAX MIN ) q 

[16] *0 a 

[17] p 


NORMALIZE REGION COORDS. 
RESHAPE NORMALIZED COORDS. 
NORMALIZE XY WITH MAX MIN 
RESHAPE NORM. RECION X Y 
FIND CLOSEST REGION 
DO REGRESSION ANALYSIS 
UNNORMALIZED THE RESULT 
EXIT 


[18] NOXY:R*-(( 2 + 1 + pD ) , 3 , 0 )p 0 p 


RETURN IP NO VALUES 


INTQRLREG 


[ 0 ] 

[ 1 ] 

[ 2 ] 

[3] 

[4] 

[5] 


R+D INTQRLREG XY ; 4 ; B ; E : I ;NDV ;RI ;UI 

p CALLED BY INTQR, IT PERFORMS DOUBLE QUADRATIC REGRESSION BY RECION 
p XY IS A 3 COLUMN MATRIX OF X AND Y VALUES AND REGION NUMBERS 
p D IS RECION X Y Z DATAp (S3) ( NUMBER OF REGIONS ) 8 
* R IS INTERP. VALUES AND 2 PART. DERIV. p r2+ - l tp = D [ 2 ] ) 3 (ltpXY) 


[6] RI-*XY [ ; 3 ] p 

[7] XY«-XY[;1 2] p 

[8] E+3 1 2 <S)D a 

[9] UI-*REMDUPEL RI a 

[10] NDV-*(~2 + ~l*pE) p 

[11] R* (NDV , 3,l+pXY)pO p 

[12] A*- (NDV >P 4)p4«-((2 + p£ , ) f 2)+ff p 

[13] E*-A , [4] 2 3 l$((2tp£),-tfCJO+£ P 

[14] p 

[15] p 


RECION INDICIES 
X AND Y DATA 
TRANSPOSE REGION DATA 
UNIQUE RECION INDICIES 
NUMBER OF DEPENDENT VARIABLES 
INITIALIZE RETURN VARIABLE 
RESHAPE INTERP. DATA ( INDEP . VARS. 
ADD BACK ACAIN THE DEPEND. VARS. 

E HAS p NDV ( NBP. . RECIONS ) 8 3 


) 


[16] LOOP:-*- (0=pUI )/0 p 

[17] 4+c [2 3 ] £T [ ; U I [ 1 ] ; ; ] q 

[18] I+IEAl UIZ12-RI a 

[19] B+c [ 2 3]((p4),pB )pB-*-XY[ J : ] p 

[20] 4«o4 QUADRGXY-B p 

[21] B*MINMAXt)2 0 *D [;£//[ 1] ; ] p 

[ 22 ] 4[; ;1]*H(I»4[; ;1] )LIMITT0 B a 

[23] fl[;;I]t-l 3 2$4 p 

[ 24 ] UI+1 + UI a 

[25] +L00P p 


LOOP ON EACH REGION OF INTEREST 
NESTED ARRAY OF RECION DATA 
INDICIES OF XY WITHIN THIS RECION 
NESTED ARRAY OF XY IN THIS RECION 
DBL. QUAD. REGRES. ON “ DEP . VAR. 

MIN AND MAX VALS. FOR REGION p 2 NDV 
LIMIT INTERP. Z VALS. TO REGION Z'S 
ASSIGN RETURN VALUES 
DROP ONE FROM LOOP VARIABLE 
END OF LOOP 
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INTQRLUNN 
CO] 

Cl] 

C2] 

C3] 

C4] 

C5] 

C 6 1 
C7] 

C8 ] 

C9] 

CIO] 

Cll] 


U+INTQRLUNN A ;R; MMM _ _ 

UNNORMALIZE A MATRIX -A- ( NESTED ARRAY.) WRT . NANCE 
AC1] NORMALIZED MATRIX BETWEEN 0 AND 1 EACH COLUMN 
A C2] MAXIMUM VALUES FOR EACH COLUMN 

AC3] MINIMUM VALUES FOR EACH COLUMN h yy1 

R RETURNS INTERP . VALUES WITH 2 PAPT. DERIV . pCDEP. 7A/?S- ) 3 (1 + pXY) 

’ / P ujj v wra/WA p MATRIX. MAX . AND MINS . 

p MIN VALUES 

MMMMPAX ( 1 ) ( C 2 +MMM ) *0PAX ( 1 )P C ! . 1 ; 3 ) « NORM . DEPENDENT VARS . 

U+U , C2] ( 2+MMM)*0PAX ( 1 ) (PC ; »2 ; ] *0PAX(2 )MMM[1] ) « y 

U+U, C2] (2+MMM)x0PAX(l)(PC; t 3;]+0PAX(2 )MMMC2] ) P VPr. X 


INV 

CO] 

Cl] 

C 2 ] 

C 3 ] 
C41 
C 5] 
[6] 

C 7 ] 
C8] 

C 9] 
CIO] 
Cll] 
C12] 


p 1 CALCULATES INVERSE OF MULTIPLE 2*2 MATRICES? 

p 


P+pA 

A-(((l+p(,A))*4),4)p,A 
I<-(pA)pl 001 

A [ ; 4 ]«- A [ ; 4 ] - A [ ; 3 ] * 4 [ ; 2 ] M [ ; 1 ] 
j[;3]<- - l*/C;l] x A[;3]*AC;l] 
j[ : i]*(JC ; 1] -II ; 3] «A [ ; 2] +AC ;4] )+A[ ; 1 
J[!23 + (IC5 2]-JC:»*]*^C5 2]*ACs 4] )+A[;l 

i’C;33*jrC«33*ACj43 

/C;4]*IC;4]+A[;4] 

I+Rf>I 


1 

1 


N 2 2 


JACCHK 

[1] p CHECK IF DETERMINANT OP JACOBIAN HAS A SIGN REVERSAL 

[2] -*-(02(L/,A)xr / , A ) /WARN 

[4] WARN: 'WARNING'. I 1 DETERMINANT OF JACOBIAN HAS A SICN REVERSAL 1 
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JACOB 

CO] JAC+RZ JACOB S;RSiDX;DE;RN-,ZN;A 

[1] a CALC. THE JACOBIAN MATRICES FOR ALL ELEMENTS 

[2] P S ARE THE SHAPE FUNCTIONS 

[3] o RZ ARE THE COORDINATES OF THE ELEMENTS p2 NE ( NODES PER ELEMENT) 

[4] p 

[5] RS+ZpS ) [2] 

C6] DX+2 1 3$(W£,~2+pS)pS[2; :] 

[7] DE+2 1 3*<tf£,"2+pS)pS[3; ;] 

C8] RN+(RS ,~2+pRZ)pRZll ; ; ] 

[9] ZN+(.RS , ~2+pRZ)pRZZ2 ; ; ] 

CIO] A+((RS*NE),l)p+/DXxRN 
Cll] A+A,((.RS*NE),l)p+/DX*ZN 

[12] A+A, ((.RS*NE) ,l)p+/DE*RN 

[13] A*-A,(.(.RS*NE),\)p+/DE*ZN 

[14] JAC+(.RS,NE,2,2)pA 


LIMITTO 

[0] R*- Y LIMITTO X 

[1] a LIMIT Y TO RANCE BETWEEN DEFINED BY X 

[2] a THE FIRST DIMENSION OF X MUST BE p2 THE MIN AND MAX VALUES 

[3] p R RETURNS MODIFIED Y p ( pY ) 

[4] p 

[5] X+MINMAX X p MAKE SURE MIN VALUES ARE FIRST 

[6] Y*YZOPAX{ppY ) ( , (1 , 1 + pX ) »X ) p SET MIN VALUES 

[7] R+Y10PAX (ppY) ( , ( l,l+pX)tX) p SET MAX VALUES 


LINFV 

[ 0 ] 

[ 1 ] 

C2] 

[3] 

[4] 

[5] 

[ 6 ] 

[7] 

[ 8 ] 

[9] 

[ 10 ] 
[ 11 ] 
[ 12 ] 


R+LINFV F-,SHP;A 

p CONVERTS LINEAR FIELD VARIABLE TO VALUES AT ALL NODES 
p F IS LINEAR FIELD VARIABLE p {ND+ 1)*2 
p R IS NODAL VALUES ( GLOBAL BASIS) 
a 


R*(NNC 2 )p0 p 

A+2 5p0 10 _ 10"10100 p 
A*-{LSHAPE A)[l;;] p 
SHP+((NE),pA)pA p 

A*-+/SHP*2 3 l*(l*pSHP)p1iF->-FZLN0DEl A 
A*-F ,A p p 

Cil 5263748 9] p 
Rl,ELN)*,A p 


INITIALIZE RETURN VARIABLE 
MID NODES XI ETA COORDS. 

SHAPE FUNCTIONS FOR MID NODES 
RESHAPE FOR ALL ELEMENTS 
MID NODE VALUES 
COMBINE WITH LIN. NODE VALUES 
REORDER ELEMENT BASIS 
RETURN VARIABLE 
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LSHAPE 

[0] SHP+LSHAPE XE;C;ETA;XI 

[1] a CALCULATES LINEAR SHAPE FUNCTIONS 

[2] A 

[3] XI«-X£[1:] A 

[4] ETA*-XE12 ;] a 

[5] C«-l ,ETA ,XI , Cl . 51ETA*XI A 

[6] C* 1 3 2 4*(3,4,pC)pC a 

[7] SHP++/C* 2 1 3 4*((ltpX.T),p£££)p£££ 


AND THEIR DERIVATIVES 

XI COORDINATES 
ETA COORDINATES 
POLYNOMIAL COEFFICIENTS 
RESHAPE TO MATCH Q$L 
a CALCULATE SHAPE FUNCTIONS 


LUMP 

[0] R*LUMP C 

[1] a LUMPS THE CAPACITANCE MATRIX 

[2] A 

[3] C++/C A 

[4] R*-DIAC C a 


BY ROWWISE SUMMATION 

ROWWISE SUMMATION 

USE SUM TERMS FOR DIAGONALS 


LX 

[0] LX ; A 

[1] a LATENT EXPRESSION 
[23 A 

[3] STACK *- 1 CMS ( 1 9 2 1 a 

[4] A«-10 0 □ SVO ' ' a 

[5] 4 -*-101 □ SVO ' ' « 

[6] STAQK*-' )EDITOR 2' a 

[7] CLS a 

[8] STRIPE a 

[9] SKIP 1 a 

[10] ' PHASTRAN ' a 

[11] SKIP 1 a 

[12] W£ID*-' PHASTRAN' a 


INITIALIZE STACK VARIABLE 
SHARE SYSTEM VARIABLE 
SHARE STACK VARIABLE 
USE EDITOR 2 
CLEAR THE SCREEN 
DISPLAY A STRIPE ON SCREEN 
SKIP A LINE 
HEADER 
SKIP A LINE 
WORKSPACE NAME 


MAP 

[ 0 ] 

[ 1 ] 

[ 2 ] 

[3] 

[4] 

[5] 

[ 6 ] 

[7] 

[ 8 ] 

[9] 

[ 10 ] 
[ 11 ] 


MAP; Cl ETA ; X ; SHP ;XI;Y 

A MAP XI-ETA COORDINATES OF NODES INTO X-Y SYSTEM 


XI+XENi 1;] A 
ETA*-XENL 2;] a 

C*(ETAxXI*2 ) , (XI*ETA*2 ) , [1 . 5] (XI*2 )xETA*2 A 
C*\,XI,ETA , (XI » ETA ) , (XI* 2 ) , (ETA*2 ) ,C a 
SHP*-+/CxOPAX( i 3 )<((i+pC),p£Sfl££>p£££££) « 
X*-+/SHPx0PAX(2)RZRLUl a 
Y++/SHPxOPAX(2)RZR(2;l a 
flZ/V«-X,[0.5]Y a 
RZE+XLELN1 ,[0.5]Y[E£tf] a 


XI COORDINATES 
ETA COORDINATES 
BUILD MATRIX COLUMNS OF 
.XI-ETA POLYNOMIAL 
SHAPE FUNCTIONS 
X COORDINATES 
Y COORDINATES 
COORDS. OF GLOBAL NODES 
COORDS. OP ELEMENT NODES 
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MDET 

[0] A+MDET B 

Cl] o CALCULATES DETERMINANT OF MULTIPLE 2*2 MATRICESd N 22 

[2] R 

[3] A *( ( ( 1 tpA )*4 ) , 4 )pA* ,B 

[4] A«-UCjl]*A[:4])-iiC;2]*A[s3] 

[5] d«-(2+((ppfl)-2)4*pfl)pA 


MEAN 

[0] M+MEAN A 

Cl] R FINDS MEAN OVER LAST INDEX OP A VECTOR 
C2] M*- ( *~1 + pj4 )x + /A 


M INMAX 

CO] R+MINMAX X 

Cl] R RETURNS THE MINIMUM AND MAXIMUM VALUES OF X 
C 2 ] fl«-(L/X),[0.5]r/X 


NNC 

CO] N--NNC FT \NLN \NQN 

Cl] R RETURNS THE TOTAL NUMBER OF NODES (.GLOBAL) FOR EACH FIELD TYPE 
C2] r 

C3] NLN+(ND+1)*2 a NUMBER OF LINEAR NODES 

C4] NQN+ (1 + 2 *ND ) * 2 fl NUMBER OP QUADRATIC NODES 

C5] N+(.NLN ,NQN)IFELTIFT11 a RETURN VARIABLE 
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NODEM 
CO] 
Cl] 
C2] 


R+NODEM A : B1 :B2 ; C iDiRtRNODE i CNODE 

assemble hooel matricies from element matricies 

a 


[3] 

S+ltltpA a 

[4] 

O" 1 1 1 + P* A 

[5] 

+(S=4)/SI A 

[6] 

RN0DE*ELN a 

[7] 

B\*NNG 2 a 

[8] 

+CTYPE A 

[91 

RL : RNODE+LNODE a 

[10] 

Bl+(.ND+l )*2 A 

[11] 

CTYPE:+(C=*)/CL a 

[12] 

CNODE+ELN a 

[13] 

B2*NNG 2 a 

[14] 

+NXT a 

[15] 

CL : CN0DE*-LN0DE a 

[16] 

B2*-(.ND + \ )*2 a 

[17] 

NXTlD+2 3 l*l<.R,C,NE)pS)CNODE a 

[18] 

D«-,C + 3 2 l$(C ,R ,NE)p$(RN0DE- 

[19] 

D*D REDUCE, A a 

[20] 

fl«- ( x /B1 , B2 )p 0 a 

[21] 

BCDC2;]]«-,0Cls] a 

[22] 

B-<-(Bl ,B2 )pB a 

NODEV 

[0] 

B-NODEV A;B1 ;D;R\ RHODE 

[1] 

p ASSEMBLE NODEL VECTOR FROM 

[2] 

0 

[3] 

S«-l + 1 + p A a 

[4] 

-*-(S = 4)/SL A 

[5] 

RNODE+ELN a 

[6] 

B\*-NNG 2 a 

[7] 

->-NXT a 

[8] 

RL : RN0DE*-LN0DE a 

[9] 

BH-CSD + 1 )*2 a 

[10] 

NXT'.D* ,RN0DE a 

[11] 

D*-D REDUCE, A a 

[12] 

B«-Blp0 a 

[13] 

BC0C2;]]-*-DCl;] a 

[14] 

B*(Bl,l)pB a 


NUMBER OF ROWS 
NUMBER OF COLUMNS 
CHECK IF LINEAR ROWS 
USE QUADRATIC NUMBERING 
NUMBER OF QUADRATIC NODES 
JUMP AND CHECK COLUMN TIPS 
USE LINEAR NUMBERING 
NUMBER OF LINEAR NODES 
CHECK IF QUADRATIC COLUMNS 
USE QUADRATIC NUMBERING 
NUMBER OF QUADRATIC NODES 
JUMP 

USE LINEAR NUMBERING 
NUMBER OF LINEAR NODES 
RESHAPE COLUMN NODE NUMBERS 
COMBINE WITH ROW NUMBERS 
SUM WHERE MULTIPLE INDICIES 
INITIALIZE RETURN VARIABLE 
REORDER RETURN VARIABLE 
RESHAPE RETURN VARIABLE 


NUMBER OF ROWS 
CHECK IF LINEAR 
QUADRATIC NUMBERING 
NUMBER OF QUADRATIC NODES 
JUMP 

LINEAR NODE NUMBERING 
NUMBER OF LINEAR NODES 
STRING OUT NUMBERS 
SUM WHERE MULTIPLE INDICIES 
INITIALIZE RETURN VARIABLE 
REORDER RETURN VARIABLE 
RESHAPE RETURN VARIABLE 


NORMR 

CO] 

Cl] 

C2] 

C3] 

C4] 

C5] 

C6] 

C7] 

C8] 

C9] 


'.+NORMR FiMIN-.MAX ,„„„ rir „ ro . 

NORMALIZE MATRIX ”F" WRT . RANGE ( RESULT IS A 
SCI] NORMALIZED MATRIX BETWEEN 0 AND 1 
SC2] MAXIMUM VALUES FOR EACH COLUMN 
SC3] MINIMUM VALUES FOR EACH COLUMN 


NESTED ARRAY ) 


MAX+t/F « 

MIN+ISF a 

R*-(F- (f>F)pMIN)* (pF)pMAX-MIN a 
R+R MAX MIN a 


IXIMUM OF EACH COLUMN 
NIMUM OF EACH COLUMN 
)RMALIZE FROM 0-1 
■.TURN NEW ARRAY , WITH MAX. 


AND MIN 
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NS 

[ 0 ] 

Cl] 

[ 2 ] 

[3] 

[4] 


R+NS T 

e RETURNS SHAPE FUNCTIONS FOR ELEMENT TYPES 
« T IS THE FIELD TYPE ( E.C . 1-ENERGY ) 
a 

R+^NSTlFELTim a SELECT FROM GLOBAL NESTED ARRAY 


OPAX 

CO] R+X (FN OPAX )B ; AS ;MIL; RT;Y 

Cl] a APPLIES DIATIC PRIMATIVE OPERATOR ALONG AXES 


1 2 j a 

C3] (AS Y)*B a 

C4] *((ppX)=ppY)/APP a 

C5] +( (ppY)>ppX)/YL a 

C6] MIL+~(\ppX)eAS a 

C7] *(-^/(pY)=(-MIL)/pX)/ERRl a 

C8] RT--(MIL/\ppX) , (*MIL)/\ppX a 

C9] y-flr«l((pX)CAr] )py a 

CIO] +APP a 

Cll] YL:MIL*~(\ppY)eAS a 

C12] + (~*/(pX')=(~MIL)/pY)/ERRl a 

C 1 3 ] RT*(MIL/\ ppY ) , (-MIL )/ i ppY a 

C 1 4 ] X*RT<*((pY)lRT'\ )pX a 

Cl 5] APP-.R+X FN Y a 

C16] +0 a 

Cl 7] ERR1-.DES 'RANK ERROR' 


RIGHT ARGUMENT CONTAINS AXES AND Y 

IF RANKS ARE SAME APPLY THE OPERATOR 

FIND LARGEST RANK 

MISSING INDEX LOCATIONS 

CHECK FOR RANK ERROR 

RESHAPE AND TRANSPOSE INFO . 

RESHAPE Y 

JUMP TO APPLY THE OPERATOR 
MISSING INDEX LOCATIONS 
CHECK FOR RANK ERROR 
RESHAPE AND TRANSPOSE INFO. 

RESHAPE X 

PERFORM THE OPERATION 
EXIT 


PH AST RAN 
CO] PHASTRAN A 

Cl] a MAIN CONTROL FUNCTION FOR PHASE 
C2] a 

C 3] FRONTEND A a 

C 4 ] TIME STEP a a 

C 5] RESULTS* RESULTS a 

C6] STATUS TSTOP a 


CHANGE ANALYSIS MODEL 

UPFRONT, ONCE ONLY FUNCTIONS 
TIME INCREMENT FUNCTION 
FORMAT FINAL RESULTS 
DISPLAY COMPUTER USAGE 


PROPDA TALSEL 


CO] 

Cl] 

C2] 

C3] 


R+PROPDATAbSEL N 

a RETURNS SUBSET OF PRQPDATA BASED ON STATE CONDITIONS AND PROPERTIES 
« N IS NUMERIC VECTOR CORRESPONDING TO SF NOTATION 


C4] R*PRQEQAT& a 

C5] fl+(v//?C3: j]eJV)/C2]/? a 


ALL PROPERTY DATA 
ONLY REQUESTED STATES 
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PRSCRB 

CO] 

[ 1 ] 

[ 2 ] 

[3] 

C4] 

[5] 

[ 6 ] 

[7] 

[ 8 ] 

[9] 

CIO] 

[113 

[ 12 ] 

[13] 

[14] 

[15] 


PR $ CRB * A * B NP * K.S 

n prescribes finite element equations in The global vars. kbar 

p PP£Q[1;3 POSITIONS IN KBAR (1 THRU 1 + p KBAR) 

p PFEQ [ 2 ; ] CORRESPONDING PRESCRIBED VALUES 

p PFEQ13 : ] (0 FOR SOLID NODES, 1 FOR BOUNDARY CONDITIONS ) 


AND RBAR 


PFEQ+PPEQ [ i &PFEQ [ 1 ; ] ] p 
PFEQ*-(RDUPL PFEQLi.il )/PFEQ p 
A+\l+pKBAR p 

PFEQ *■ ( PFEQ [ 1 : ] e A ) /PFEQ p 
BNP*-~AfPFEQ C 1 ;] p 
A+-(~BNP)/KBAR a 
KB *- ( ( (>A)pPFEQL2ll )*A « 

KB<- + /Pf’£'G[3;]x[2]KB n 
RBAR*-(BNP/RBAR)+BNP/KB n 
KBAR+8NP/ 111 BNP/ KBAR o 


PUT IN ASCENDING ORDER 
REMOVE DUPLS. ( LEAVE 1ST) 
ALL POSITIONS IN KBAR 
REMOVE ANY OUT OF RANGE 
BOOL. POS. NOT PRSCRBD. 
COLUMNS OF K'S PRESCRIBED 
K'S * PRESCRIBED VALUES 
ZERO SOLID NODES 
NEW R VECTOR 
NEW K MATRIX 


PRSIDE 

[0] ' PRSIDE FT -.B-.C-.D 1;D2 

[1] o MODIFIES PFEQ WHICH DEFINES 

[2] p 

[3] B*BCl; 3 2 1 4 5] p 

[4 3 B*-RN FSTCM B p 

[5] fl«-l FSTCM B a 

[63 ( (pFDpcfl )PRSIDELFLD"FT a 

[73 ENDxDl+PFEQZlxl p 

[83 Dl*-<t>(C*-RDUPL Dl)/Dl*-$D1 n 

[93 D2*-$C /QPFEQ 12 3;] p 

[103 PFEQ*-D1, [1302 p 


PRSIDELADJ 

[03 AN*rON PRSIDELADJ TYPE ;PFTS 
[13 p ADJUSTS ORIGINAL NODE POSITIONS 
[23 P 

[33 PFTS*-(-l ++/*\TYPE*$FT)*FT p 
[43 AN+ON++/NNC PFTS p 


PRESCRIBED BOUNDARY CONDTIONS 

REORDER BC (.REG. TYPE FT SIDE VAL.) 
ONLY THIS RECION 
ONLY PRESCRIBED BC'S 
PRESCRIBE BC'S FOR EACH FIELD TYPE 
MODIFY PFEC TO ELIMINATE 
-DUPLICATE NODE BC'S 
-LEAVING ONLY THE FIRST 
RESET GLOBAL VARIABLE PFEQ 


IN A MATRIX FOR FIELD TYPE 

PREVIOUS FIELD TYPES 
ADJUSTED NUMBERS ( POSITIONS ) 


PRSIDELFLD 

[03 B PRSIDELFLD TYPE;N-,V 

[13 P PRESCRIBES BOUNDARY CONDITIONS 

[23 p 

[33 B+TYPE FSTCM B R 

[43 -*-(0 = pB )/0 p 

[53 N*TYPE SIDENODES B[:13 p 

[63 y«-,$(<l>pA0p«"»"B[ :23 p 

[73 (N V )<-((, N)V)PRSIDELTEMP TYPE P 

[83 N*-(.N)PRSIDELADJ TYPE p 

[93 PFEQ*-PFEQ,N, [13 V, [0. 531 p 


FOR A FIELD 

BOUNDARY CONDITIONS FOR TYPE 
EXIT IF NO BC'S FOR THIS TYPE 
FIND NODE NUMBERS 
FIND ASSOCIATED VALUES 
HANDLE TEMP. BC'S DIFFERENTLY 
ADJUST MATRIX POSITIONING 
ADD PREVIOUS PRESCRIBED VALUES 
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PRSIDELTEMP 

CO] D+NV PRSIDELTEMP TYPE;I:CiIl:I2:N:V 


Cl] o CONVERTS TEMPERATURE BOUNDARY 

[2] a 

[3] D*NV a 

[4] +(TYPE*1)/ 0 a 

C5] (« V)*NV a 

C6] 0+A7, [0. 5] (piV)pO o 

[7] n*IEAl I+( 1 FVCB SP)CW]el 4 6 

[8] I+lPROPDATALSEL 1 4 6) [4 1 2 ; ; ! 

[9] C+/CJ1] ,C1.5](1 PVCB PF)iNUW. 

CIO] 0C2;/1]«-(C INTQR 1)1 1:1;] a 
Cll] D*-c [ 2] 0 o 


CONDITIONS TO ENERGY BC'S 

INITIALIZE RETURN VARIABLE 
EXIT IP NOT TEMPERATURE BC 
NODES AND VALUES 
SETUP RETURN VARIABLE 
o IND. OP NODES IN SINGLE PHASE 

« SINGLE-PHASE PROPERTY DATA 

« USE TEMPERATURE AND PRESSURE 

INTERPOLATE SINCLE-PHASE POINTS 
RETURN NESTED ARRAY 


PRSPRES 
CO] PRSPRES 

Cl] a PRESCRIBES PRESSURE NODEM 
C 2] a 

C3] -*■ ( 2 = DNC 'PRNODE' )/NP a 

C 4 ] 0 

C 5] PRSP 1 a 

C6] -»0 

C7] NP-.PRSP 2 a 


CHECK IF PRESSURE NODE PRESCRIBED 
PRESCRIBE PRES. ( TOTAL MASS CONSTRAINED) 
PRESCRIBED PRESSURE NODE ( OPEN SYSTEM) 


PRSP 1 

CO] PRSP1 ;A;B;NN 

Cl] a PRESCRIBES PRESSURE NODEM : CASE 
C2] a USES PRESSURE VALUE PROM LAST 
C3] a 

C4] NN+NNG 2 a 

C5] A+NNp ( ( 1 +2 *ND )pl 0),(WD+l)p0 a 
C6] B + v /(2 ,NN)p(.NN+\2xNN)ePFEQLl;) a 

C 7] /J«-U*0 (ND + 1 )*2 a 

C8] B*PP[U-A) a 

C9] PFEQ*PFEQ, 3 lp ( ( 1 1<4 ) + 3*AW ) ,S , 1 a 


WHERE TOTAL MASS CONSTRAINED 
ITERATION 

NUMBER OF QUAD. NODES ( GLOBAL BASIS) 
BOOL. WHERE PRES. NODES OCCUR 
BOOL. WHERE VELOCITIES PRESCRIBED 
PRES. NODES WHERE VEL . NOT PRSCRBED . 
USE 1ST UNPRSCRBED . PRES. NODE 
PRESCRIBE THE PRESSURE NODE 


PRSP 2 

CO] PRSP2 ; N 

Cl] a PRESCRIBES PRESSURE NODEM: CASE OP OPEN SYSTEM 
C2] a 

C3] N+(+/NNG 2 3 )+PRN0DEll] a MATRIX LOCATION OF NODE NUMBER 

C4] PFEQ+PPEQ , 3 IpN .PRNODEL2 ) ,1 a PRESCRIBE THE PRESSURE NODE 
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PRSVLC 

[0] CE+PRSVLC FT 

[1] a RETURNS LOCATIONS OF 
C2] A 

[3] CE+ (+/NNG FT ) pi a 

[4] CEUPFBQZlil )3<-0 a 


PRESCRIBED VALUES FOR THESE FIELD TYPES 

VECTOR OP l'S FOR ALL NODES 
INSERT 0'S WHERE PRESCRIBED 


QLCSRP 

[0] SHP+QLCSHP XE;A;C-.ETA;XI 

[1] a CALCS. QUADRATIC LAGRANCIAN SHAPE FUNCTIONS 

[2] R 

[3] XI-<-mi;] R 

[4] EZM+XEC2;] R 

[5] C+(ETA*XI*2 ) , (XI*ETA+2 ) , Cl . 5] (XI* ETA )*2 R 

C6] C+l.XI ,ETA,(XI*ETA) ,(XI*2) ,(ETA*2) ,C a 

C7] SHP+C+.xQCSQLQ R 

C83 A+QSQLQxOPAX 2(010120212)r 

C9] i5-*-i4C;2 5471839 6] R 

CIO] SHP+SHP ,C0.5]C+.*$4 n 

Cll] A+C£QLG*OPAX 2(001102122)r 

C12] ,4-w5C:3 4687192 5] R 

C13] SHP+SHP, [1]C+.»W R 


AND DERIVATIVES 

XI COORDINATES 

ETA COORDINATES 

BUILD MATRIX OF POLYNOMIAL 

. COEFFICIENTS 

CALCULATE SHAPE FUNCTIONS 

COEF. FOR DERIV. WRT XI 

REORDER 

DERIV. WRT XI 

COEF. FOR DERIV. WRT ETA 

REORDER 

DERIV. WRT ETA 


QUADAPLY 

CO] P+QUADLPLY XY-.X-.Y 

Cl] R FORMS POLYNOMIAL FOR 2 DIMENSIONAL QUADRATIC 
C 2 ] fl XX IS A 2 COLUMN MATRIX OF X Y VALUES 
C 3 ] R 

C4] x<-xxC:i] r 

C 5] Y«-XYC ; 2] r 

C6] P+1 ,X, (X*2 ) ,Y , (Y*2 ) , (Xxy ) , (YxX*2 ) , Cl . 5] (XxY*2 ) 


X VALUES 
Y VALUES 
FORM POLYNOMIAL 


QUADRCXY 


CO] 

Cl] 

C 2 ] 
C3] 

C4] 

C5] 

C6] 

C7] 

C 8] 
C9] 
CIO] 
Cll] 
C12] 
CIS ] 
C14] 
C 1 5 ] 


Z+D QUADRCXY XY ; C ; K ; DX ; DY : PXY 

a PERFORMS QUADRATIC REGRESSION ANALYSIS IN 3 DIMENSIONS X Y Z 
a XY IS A 2 COLUMN MATRIX OF X AND Y VALUES TO BE INTERPOLATED 

a D IS A 3 COLUMN MATRIX OF X Y Z DATA USED IN THE INTERPOLATION 

a Z RETURNS INTERP. VALUES AND THEIR DERIVATIVES p (lfpXX) 3 


K+QUADLPLY D Csl 2] R 
K+DZ : 3] IK R 
PXY+QUADLPLY XY fl 
Z-PXY+.xK R 

C+ 1 201120 0xKC2 3168711] a 
DX+PXY + . *C R 
Z+Z , Cl . 5]DX a 

C+ 1 112020 0xXC4 6751811] R 
DY+PXY+.xC A 
Z+Z.DY a 


POLYNOMIAL OF INTERP. DATA 
REGRESSION ON DATA 
POLY. OF XY'S TO BE INTERP. 
INTERPOLATED VALUES 
COEF. FOR DERIV. WRT. X 
DERIV. OP Z WRT. X AT XY 
CATENATE TO RETURN VARIABLE 
COEF. FOR DERIV. WRT. Y 
DERIV. OF Z WRT. Y AT XY 
CATENATE TO RETURN VARIABLE 
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RCOND 


[0] 

RC*RC0ND 



Cl] 

fl RETURNS SUB -VECTOR FOR INTERNAL CONDUCTION 



[23 

fl 



[3] 

flO+/[3](OX 1 ),DY 1 )*0PAX ( 2 3)(1 FVEB TF) a 

DERIVATIVES OF TEMP. 

[4] 

RC*RC*0PAX{\ 2 ) {FVIP XT) a 

MULT. BY 

CONDUCTIVITY 

C 5] 

RC**/RC*0PAX { 1 2 4 ) ( OX 1),DY 1) n 

MULT. BY 

DERIV. SHP. FNS. 

[6] 

RC*{{pRC) ,l)pRC R 

RESHAPE 

TO COLUMN VECTOR 

[7] 

RC*-N0DEV AINTCRT RC fl 

I NT CRT. 

GLOBAL NODE BASIS 

RDUPL 



C03 

R+RDUPL X 



Cl] 

R RETURNS A BOOLEAN FOR REDUCINC DUPLICATE VALUES LEAVING ONLY THE 1ST 

C 2 ] 

a 



C 3 ] 

R + 1 V*<\X°.-X 




REDUCE 


[03 

X*IND REDUCE VAL ; C ; D ;C2 ;C3 ; I 


C13 

fl REDUCES {BY SUMMATION) A VECTOR 

WITH MULTIPLE INDICIES 

[2] 

R 


C 3 ] 

IND*INDU*kIND'\ a 

REORDER INDICIES SEQUENTIALLY 

C4] 

VAL+VALII1 r 

AND THE CORRESPONDING VALUES 

C 5 ] 

C2*INDxl<S>IND fl 

LOCATIONS OF 'DUPLICATES 

C6] 

-*■( 1**/C2 )/NXT fl 

JUMP IF DUPLICATE INDICIES EXIST 

C7] 

X*VAL,tO .b)IND a 

RETURN VARIABLE 

[8] 

+ 0 R 

EXIT 

C 9 ] 

NXT: 


[10] 

C3*C2*(~C2 )* ( (pC2 )-l fl 

BOOLEAN WITH 0 ' S AT FIRST DUPL. 

Cll] 

D*C3 / INDxC3 fl 

REDUCED SET OF INDICIES 

[12] 

C*-VAL + l<f>VAL*(-l+pVAL)‘t>C3xl fl 

SUM VALUES 

[13] 

C*({~l*pVAL)<t>C3)/C fl 

ELIMINATE THOSE VALUES 

[14] 

X*D REDUCE C a 

RECURSION TO FURTHER REDUCE 

REMDUPEL 


CO] 

R*REMDUPEL X 


Cl] 

fl REMOVES DUPLICATE ELEMENTS OP X 


C 2 ] 

R 


C 3 ] 

R* ((XiX)=ipX)/X 


RESULTS 


CO] 

R*RESULTS ; HDR 


Cl] 

fl DISPLAYS FIELD VARIABLES 


[2] 

R 


C3 ] 

R*PF , EF , TF , SF , UF . [ 1 . 5 ] VF 



[>♦] HDR*' PRESSURE' 'ENERGY' 'TEMPERATURE' 'STATE' ' U-VELOCITY ' 'V-VELOCITY' 
C5] R+9HDR.111R 
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RBOZERO 

[0] R*-RH0ZER0;DATA’,IP 

[1] a CALCULATES THE REFERENCE DENSITY 

[2] A 

[3] DATA*(.PROPDATALSEL 4)[1 4 5;;] a 
C4] IP*- 1 2p (MEAN PF),MEAN TF A 

[S3 R*-1*,IP INTQR DATA A 


USE LIQUID DATA 

USE MEAN PRES. AND TEMP. 

INTERPOLATE ON DENSITY 


RPRES 

[0] RP*-RPRES 

[1] a FORMS PRESSURE RESULTANT VECTOR 

[2] A 

[3] RP*-PFZELN1 « (DFN UF) [ 1 ; : ] + (DFN VF) C2::3 
Cu] RP*--N0DEV AINTGRT (NS 1)*C2 3 1RP 


SAVEFIELDS 

[03 NTS SAVEFIELDS SFN „„ 

[1] a SAVES FIELDS AT EVERY DFN TIME STEP, PUTS IN NESTED 

[23 a 

[33 +(0*SFN\NTS)/0 a 

[43 £AVEFLD£*-£AVEFLD£,zTIME(PF ,EF ,TF ,SF ,UF , [1.53VF ) a 


ARRAY SAVEFLDS 

EXIT IF NOT 
SAVE FIELDS 


SET LA ST 
[03 SET LAST 

[13 a SETS VARIABLES FROM LAST 
[23 A 

[33 RELAST*-EREQ a 

[43 RFLAST*-FREQ a 

[53 EPHI*-EF a 

[63 FPHI*-UF ,VF , 4 FVCB PP a 


SETSOLID 

[03 SETSOLID \A\N\ NSOL j P 
[13 a SETS VELOCITIES TO ZERO 
[23 a 

[33 NSOL+IEAl~SFe ' + a 

[43 UFtNSOLl+O a 

[53 VFZNSOLl+O a 


ITERATION 

LAST ENERCY RESULTANT 
LAST PLOW RESULTANT 
LAST ENERGY SOLUTION 
LAST FLOW SOLUTION 


NODES IN SOLID STATE 

FIND NON-LIQUID NODES 
SET U VELOCITY TO ZERO 
SET V VELOCITY TO ZERO 
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SETSTAR 

[0] SETSTAR 

Cl] p SETS EE, UP, VP AND PF AT CAST 

[2] p 

[3] EFSTAR*EF a 

[4] +(FC= 0)/0 a 

C 53 UPSTAR*-UP a 

C 6 3 VFSTAR+VP p 

[7] PFSTAR+PF a 


CONVERGENCE ITERATION 

LAST ENERGY VALUES 

CHECK IP PLOW CALCS. WERE SOLVED 

LAST 0 VELOCITIES 

LAST V VELOCITIES 

LAST PRESSURES 


SIDENODES 

[0] R*-FT SIDENODES SiET 

[ID p RETURNS THE REGION SIDE NODES FOR A FIELD TYPE 
[2] a S IS THE SIDE NUMBERS 


[33 p FT IS THE TYPE OP FIELD ( E.G . 

[4] a 

[5] ET+FELTIFT1 a 

[6] -*-(£T=l 2)/Tl,T2 a 

[73 +0 n 

[8] Tl:R*-RSLNlS;) a 

[9] ♦ 0 o 

[10] T2 :R+RSNZSt ] o 


SI FAC 

[0] SIF+SIFAC;DX;DE;XN;YN-,A;B;C:D-,Sl\ 

[1] n CALCULATES THE SIDE INTEGRATION 

[2] o 

[3] DX *- 2 1 3S>(.NE,pQSSZ2i ;] )p«SS[2; j] 

[4] DE->-2 1 3$(.NE ,pQSSL3 ; ; ] )pQSS[3 ; ; ] 

[5] XN+( (pQSS ) [2] ,~2fpRZE)pRZEil ; ; ] c 

[6] YN+((pQSS)Z21 ,~2*pRZE)pRZEZ2; ;] c 

[7] i9-*-(4,((pflSS)[2]t4),WE)p+/ZJX*X« o 

[8] E«-(4, ((pCSS)[2] + 4),tfE)p + /i>XxyjV n 

[9] 0(4,((pQSS)[2]+4),N£:)p + /D£xXtf p 

[10] D*(4, ((pGSS)[2]+4),Wff)p+/CFxy« A 

[11] SH-,*(U[l::]*2)+fl[l;;]*2)*0.5 A 

C12] S2-,$((C[2; ;]*2)+£>[2; ;]*2)*0.5 a 

[13] 53*-,ty((4[3; i ] *2 )+B [ 3 ; s ] *2 )*0 . 5 p 

[14] S4«-,*((C[4; ;]*2)+i>[4; ;]*2)*0.5 P 

[15] A* 2 1 3*(/Vr,p<JSS[l;;])pQSS[l;;] p 

[16] A+(*,(SNH-1) ,NE)p+/A a 

[17] SIF+Ax (p/4 )pHSl ,S2 ,S3 , [1 . 5]S4 p 


i -ENERGY) 

CHANGE FIELD TYPE TO ELEMENT TYPE 
CHECK TYPE OF ELEMENT 
EXIT IF NOT VALID ELEMENT TYPE 
LINEAR SIDE NODES 

EXIT 

QUADRATIC SIDE NODES 


S2;S3 ;S4 

FACTOR FOR USE WITH SI NT CRT 

p DERIV. SHP. FNS. WRT. XI 

p DERIV. SHP. FNS. WRT. ETA 

X COORDS. OF NODES 
Y COORDS. OF NODES 
DERIVATIVES * COORDINATES 
DERIVATIVES x COORDINATES 
DERIVATIVES * COORDINATES 
DERIVATIVES * COORDINATES 
SIDE 1 
SIDE 2 
SIDE 3 
SIDE 4 

INCLUDE THE SHAPE FUNCTION 
.AS PART OP THE . . . 

SIDE INTEGRATION FACTOR 
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SINTCON 

[0] SINTCON N-,A\N 1 

[1] a INTEGRATION CONSTANTS FOR USE WITH 

[2] A 

[3] Nl+l+SNI+N a 

[u] A+N1*(QIL1; ;] KSWj] a 

[5] SXEI*A,(Nlpl) ,A,(Nlp~l) A 

[6] SXEI+SXEI , [0 . 5] (Allp“l ) ,A , (.Nipl ) ,A A 

[7] QSS+QLCSHP SXEI a 

[8] LSS+LSBAPE SXEI A 

[9] QSIF+SIFAC a 


SINTCRT 

ORDER OF GAUSS-LECENDRE 
INTEGRATION POINT CONSTANTS 
XI COORDINATES 
ADD ETA COORDINATES 
QUADRATIC SHAPE FUNCTIONS 
LINEAR SHAPE FUNCTIONS 
SIDE INTEGRATION FACTOR 


SINTCRT 

[0] Ii-SINTGRT A;N1;WE 

[1] a INTEGRATES FUNCTION OVER SIDES 

[2] A 

[3] N1+SNI+ 1 a 

[4] -»(lep .A'i/SCALAR a 

[5] -*■ ( a/ (pQSIF) = 3 tp<4 ) /NXT a 

[6] STATUS ERRMSQZ21 a 


OP ELEMENT IN XI -ETA COOR . SYSTEM 

1 +0RDER OF INTEGRATION 
CHECK IF A IS SCALAR 
CHECK IF p A IS LIKE pQAIF 
MESSAGE TO USER 

EXIT 

MULTIPLY BY INTCRT . FACTOR 
JUMP TO FINISH INTEGRATION 
MULTIPLY BY INTCRT. FACTOR 
'WEIGHTING FACTORS 
INTEGRATE ALL SIDES 


[7] +0 A 

[8] NXT:A*-Ax (.~lQ\ppA)$(~l , $pA)pQSIP a 

[9] +END a 

[10] £CALAR’.A+A*QSIP a 

[11] END:WE*-{2<t>\ ppd )$(24>Pi4 )pM+ (CI[2 : ; ] ) ISNI; ] a 

[12] ~I**/Lll*/L11WE»A a 


SKIP 

[0] S*-SKIP N 

[1] a PRODUCES "N" BLANK ROWS 

[2] S<-(AM)p' ' 


STATPRES 

[0] SPF+STATPRES-.A-.NXN-.XYC 

[1] a CALCULATES THE STATIC PRESSURE DISTRIBUTION 

[2] A 

[3] NXN+ (.NS 4 )INPROD 1 2 4 36)NS 4 a 

[4] XYC+RZBZ ; !l 3 5 7] a 

[5] A**/ ZUGRZ* [1] ((pX7C)p4 FVEB RHO)*XYC A 

[6] A*NXN INPRODZ ( 3 fpNXN) ,l)pA a 

[7] A+NODEV AINTGRT A a 

[8] SPF+-.AMNODEM AINTGRT NXN A 

[9] SPF+SPF+PRN0DEZ2'} -SPFZHPRN0DE1 a 


PRODUCT OF SHAPE FUNCTIONS 
X Y COORDS. OF CORNER NODES 
DENSITY*CRA VITY * LOCATION 
MULTIPLY BY SHAPE FUNCTIONS 
INTEGRATE ON GLOBAL BASIS 
SOLVE FOR PRESSURE 
ADD REFERENCE PRESSURE 
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STATUS 

CO] STATUS HSC 

Cl] n MANAGES STATUS MESSAGES 

C2] fl 

C3] MSC**MSC fl FORMAT MESSAGE 

C4] STATMSC+STATMSC ADDROWS MSG b RECORD IN STATMSG 

C5] MSG fl D I SPLAT THE MESS ACE AT TERMINAL 


STRIPE 

CO] S*STRIPE 

Cl] fl CREATES STRIPE LINE BORDER 
C2] S*79pCU7C145] 


TIMECHK 

[0] CK*TIMECHK 

Cl] o RETURNS A 0 IF TIME LIMIT IS EXCEEDED 
C2] fl 

C3] CK * — ( 0 .001xD/17C3]-riC3] )2TIMELIM fl CHECK RUN TIME 

C 4 ] *(CK= 1)/0 q EXIT IF OK 

C 5] STATUS ' RUN TIME LIMIT EXCEEDED < o MESSACE TO USER 

C6] o TO SET S&TIMECHK 


TIMESTEP 

CO] TIMESTEP : CIC-.N IT; NTS 

Cl] fl TIME STEPPINC FUNCTION 
C2] o 

C3] NTS * 0 b 

C 4 ] CIC * lOpO o 

C 5] L00P:TIME*TIME+TCTL[3) a 

C6] NTS*NTS+1 o 

C7] STATUS 'TIME' n+TIME) a 

C8] *(TCTL C 2] < - l +TIME ) / 0 a 

C9] *{Q-CPUCHK)/D fl 

CIO] *(0=TIMECHK)/0 fl 

Cll] NIT* AGAIN 0 o 

C12] SET SOLID fl 

C 13] CIC*1*CIC,NIT b 

C14] CIC AD J STEP NIT fl 

CIS] SETLAST o 

C 16] NTS SAVE FI ELDS SFN fl 

C17] +(CHKSS CIC)/ 0 fl 

C18] *LOOP fl 


INITIALIZE TIME STEP COUNTER 
INITIALIZE CONVERGENCE ITER. COUNTER 
LOOP ON TIME 

INCREMENT TIME STEP COUNTER 

D I SPLAT TIME 

EXIT IP AT TIME LIMIT 

EXIT IP CPU TIME EXCEEDED 

EXIT IF RUN TIME EXCEEDED 

CONVERGE AND RETURN ITERATION NUMBER 

SET SOLID NODE VELOCITIES TO ZERO 

UPDATE CONVERGENCE ITERATION COUNTER 

ADJUST TIME STEP BASED ON ITER. NBR . 

SET EF. UF, VF, PF FOR LAST TIME STEP 

SAVE FIELDS AT CERTAIN TIME STEPS 

EXIT IF STEADY STATE REACHED 

CONTINUE IN TIME 
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TOMATRIX 

[5] p “CHANCES™ SCALAR, VECTOR OR ARRAY OF RANK 3 OR HIGHER TO A MATRIX 


[2] a 

[3] M+V fl 

[ 4 ] +(2 = pplO/0 B 

[5] -*(2<pp7)/Wl B 

[ 6 ] M+(.l,pV)pV+,V " 

[7] +0 b 

[8] «i:*-«-((*/( _ l+PpP)+pt r ).‘‘l + PlOpV 0 


TSTART 

[0] TSTART 

[1] b STARTS CLOCK FOR TIME CHECKING 

[2] Tl+MJ 


INITIALIZE RETURN VARIABLE 
EXIT IF ALREADY MATRIX 
CHECK IF RANK GREATER THAN 2 
CONVERT SCALAR OR VECTOR TO A MATRIX 
EXIT 

CONVERT HIGHER RANK ARRAY TO A MATRIX 


TSTOP 

[ 0 ] 

[ 1 ] 

[ 2 ] 

C3] 

[4] 

[5] 

[ 6 ] 


R+TSTOP\T2 
a DISPLAYS CPU , 


CONNECT TIMES SINCE TSTART WAS ISSUED 


A 


T2+DAI a 
SKIP 1 a 

^(«0.001xr2[2]-riC2] ), ' SEC. CPU' a 
Rt-R ADDROw'sitO . 001 *r 2 [3] -T1 [3] ) , ' SEC. CT' 


B 


CURRENT ACCOUT INFO. 
SKIP A LINE 
DISPLAY CPU USAGE 
DISPLAY CONNECT TIME 


UPPROP 

[0] UPPROP; DATA-, IP-, DP 

fil a UPDATES PROPERTIES ON GLOBAL NODE BASIS 


C 2] 

a 

[33 

DATA+PROPDATALSEL\6 a 

[4] 

IP+PF, [1.5 ]EF fl 

[53 

DP+IP INTQR DATA fl 

[63 

SF-*-DP[l ; 1 ; ] b 

[73 

TF+DP12 ; 1 : ] a 

[83 

RHO+DPl 3}1;] a 

[93 

KT+DPlH : 1 ; ] a 

[103 

vis+djpcssis] a 

[113 

UPPROPLTEMP fl 


PROPERTY DATA 
INDEPENDENT PROPERTIES 
PERFORM INTERPOLATION 
STATE 

TEMPERATURE 

DENSITY 

CONDUCTIVITY 

VISCOSITY 

PRESCRIBED TEMPERATURES 
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UPPROPLTEMP 

CO] UPPROPLTEMP ;B;N:T 

[1] a PLACES PRESCRIBED TEMPERATURES 

C2] a 

[3] 8+BC C ; 3 2 1 4 5] p 

C4] B+RN PSTCH B a 

C5] B+l PSTCH B P 

C6] B * 1 PSTCH B a 

C7] -*(0 = ltpB)/0 n 

[8] B*»"»”fl p 

[9] N+l SIDENODES B[jl] p 

[10] r<-H($pW)pB[;2] a 

Cll] TPZ,N1+,T p 


IN UPDATED properties 

REORDER BC {REG. TYPE PT SIDE VAL . ) 

BC'S POR THIS REGION 

ONLY PRESCRIBED BC'S 

ONLY TEMP. BC'S 

EXIT IP NO TEMP. BC'S 

PRESCRIBED TEMPERATURE CONDITIONS 

BIND NODE NUMBERS 

FIND ASSOCIATED TEMP. VALUES 

REPLACE WITH PRESCRIBED TEMPERATURES 
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Cl 


267 NUMERIC 


♦+" 0.57735 

0 

" 0.33998 

0 

" 0.23862 

0 


0.57735 0 0 0 

" 0.7746 0.7746 0 0 

0.33998 " 0.86114 0.86114 0 

" 0.53847 0.53847 " 0.90618 _ 0. 90618 

0.23862 " 0.66121 0.66121 0.93247 

" 0.40585 0.40585 " 0.74153 0.74153 


0 0 
0 0 
0 0 
0 0 
0.93247 0 
0.94911 0.94911 


11 ° 
0.88889 0.55556 0 

0.65215 0.65215 0 

0.56889 0.47863 0 

0.46791 0.46791 0 

0.41796 0.38183 0 


0 0 
55556 0 0 

34785 0.34785 0 

47863 0.23693 0.23693 

36076 0.36076 0.17132 

38183 0.27971 0.27971 


0 0 
0 0 
0 0 
0 0 
0.17132 0 
0.12948 0.12948 


CSC 


4 4 NUMERIC 


+ 

++ 0.25 " 0.25 " 0.25 _ 0.25 
0.25 " 0.25 0.25 " 0.25 

0.25 0.25 0.25 0.25 

0.25 0.25 “ 0.25 " 0.25 

" 0.25 0.25 0 0 

0.25 " 0.25 0 0 

0.25 0.25 0 0 

" 0.25 " 0.25 0 0 

" 0.25 0 0.25 0 

" 0.25 0 " 0.25 0 

0.25 0 _ 0 . 2 5 0 

0.25 0 0.25 0 


Q2QLQ 9 9 NUMERIC 


0 0 0 

0 0 " 0.5 

0 0 0 

0 0.5 0 

0 0 0 

0 0 0.5 

0 0 0 

0 " 0.5 0 

10 0 


0.25 0 

0 0 
" 0.25 0 

0 0.5 

0.25 0 

0 0 
" 0.25 0 

0 0.5 

0 "1 


0 " 0.25 " 0.25 0.25 

0.5 0.5 0 " 0.5 

0 “ 0.25 0.25 _ 0 . 2 5 

0 0 " 0.5 " 0.5 

0 0.25 0.25 0.25 

0.5 " 0.5 0 " 0.5 

0 0.25 “ 0.25 0.25 

0 0 0.5 " 0.5 

"10 0 1 
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EBQEE&U 

7 3 

8 NUMERIC 





+>1000 

1000 

1000 

500 

0 

0 

0 

500 

1000 

1000 

1000 

500 

0 

0 

0 

500 

1000 

1000 

1000 

500 

0 

0 

0 

500 

2000 

3000 

4000 

4000 

4000 

3000 

2000 

2000 

4000 

4500 

5000 

5000 

5000 

4500 

4000 

4000 

5000 

6000 

7000 

7000 

7000 

6000 

5000 

5000 

1 

1 

1 

1 

1 

1 

1 

1 

2 

2 

2 

2 

2 

2 

2 

2 

4 

4 

4 

4 

4 

4 

4 

4 

“2 

“1 

0 

0 

0 

"1 

"2 

‘2 

0 

0 

0 

0 

0 

0 

0 

0 

0 

1 

2 

2 

2 

1 

0 

0 

1000 

1000 

1000 

1000 

1000 

1000 

1000 

1000 

1000 

1000 

1000 

1000 

1000 

1000 

1000 

1000 

1000 

999 

998 

998 

998 

999 

1000 

1000 

1 

1 

1 

1 

1 

1 

1 

1 

1 

1 

1 

1 

1 

1 

1 

1 

1 

1 

1 

1 

1 

1 

1 

1 

0*01 

0.01 

0.01 

0.01 

0.01 

0.01 

0.01 

0.01 

0.01 

0.005 

0.001 

0.001 

0.001 

0.005 

0,01 

0.01 

0.001 

0.001 

0.001 

0.001 

0.001 

0.001 

0.001 

0.001 
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RESULTS. 

226 6] 

PRESSURE ENERGY Tl 

0.2854 

6500 

0.1427 

6500 

0 

6500 

"0.1427 

6500 

"0.2855 

6500 

"0.4282 

6500 

*0.571 

6500 

"0.7136 

6500 

"0.8563 

6500 

"0.9989 

6500 

*1.141 

6500 

*1.284 

6500 

*1.427 

6500 

"1.569 

6500 

*1.712 

6500 

0.2854 

6012 

0.1427 

5946 

"0.0000731 

5908 

"0. 1428 

5896 

"0.2856 

5907 

"0.4284 

5945 

"0.5711 

6000 

"0.7138 

6075 

"0.8565 

6152 

"0.9991 

6238 

"1.142 

6312 

"1.284 

6376 

*1.427 

6418 

*1.57 

6439 

"1.712 

6443 

0.2854 

5632 

0.1426 

5545 

"0.0001462 

5481 

"0.1429 

5470 

"0.2857 

5502 

"0.4285 

5561 

"0.5712 

5637 

"0.7139 

5737 

"0.8566 

5852 

"0.9993 

5991 

*1.142 

6125 

"1.285 

6253 

"1.427 

6341 

*1.57 

6387 

“1.712 

6393 

0.2853 

5412 

0.1425 

5349 

"0.0002494 

5319 

"0.1431 

5339 

*0.2859 

5391 

"0.4286 

5454 

*0.5714 

5528 

*0.7141 

5612 


CHARACTER 


1.5 

1.5 

1.5 

1.5 

1.5 

1.5 

1.5 


0. 
1 . 
1. 
1 . 
1 . 
1. 


5 
5 
5 
5 
5 
5 
5 
5 

012 
0.946 
0.9075 
0.8958 
0 . 9067 
0 . 9449 
9997 
075 
152 
238 
312 
376 
1.418 
1.439 
1 . 443 
0.6318 
0.545 
0.4808 
0.4701 
0.5017 
0.5606 
0.6374 
0.7369 
0.8516 
0.9909 
1.125 
1.253 
1.341 
1.387 
1.393 
0.4121 
0.3488 
0.3193 
0.3395 
0.3905 
0.454 
0.5275 
0.612 


STATE U 
4 
4 
4 
4 
4 
4 
4 
4 
4 
4 
4 
4 
4 
4 
4 
4 
4 
4 
4 
4 
4 
4 
4 
4 
4 
4 
4 
4 
4 
4 
4 
4 
4 
4 
4 
4 
4 
4 
4 
4 
4 
4 
4 
4 
4 
4 
4 
4 
4 
4 
4 
4 
4 


V -VELOCITY 

0 
0 
0 
0 
0 
0 
0 
0 
0 
0 
0 
0 
0 
0 
0 
0 

0 . 00008336 
0 .00008165 
0.00006989 
0 . 00005228 
0.00002314 
“0 .000007881 
“0 .00003793 
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